text:
\section{Introduction}\label{sec:introduction}

A plethora of optimization-based motion planning strategies exist in the legged locomotion literature. Solving a single nonlinear trajectory optimization (TO) problem is a rather common approach. Single-optimization based methods can be categorized based on the types of dynamics representation. Reduced-order models, such as the spring loaded inverted pendulum (SLIP) model, single rigid body model (SRBM), centroidal dynamics model, etc., have been widely adopted within the community. High speed and robust running are demonstrated in \cite{di_carlo_dynamic_2018, kim_highly_2019, bledt_regularized_2020}, using SRBM and user defined heuristics. However, their framework can only generate cyclic and short-horizon gaits. To generate more complex motions, Winkler et al. proposed TOWR \cite{winkler_gait_2018}, which is capable of directly optimizing over contact sequences and timings, through its phase-based parameterization.

For reduced-order models, incorporating joint information to reason about the centroidal inertia is non-trivial, which limits its ability to generate angular-momentum-rich motions using limb movements. TOWR+ \cite{ahn_versatile_2021} improved upon the original TOWR framework by augmenting the standard SRBM with a centroidal inertia regressor network. However, it has difficulty predicting inertia trajectories for gaits not present in the training set. Wang et al. \cite{wang_unified_2021} attempted to model the leg dynamics through a point mass approximation for each limb. But, the resulting lumped-leg SRBM still cannot utilize the full capabilities of the joint motions. Earlier works, such as \cite{zordan_control_2014} and \cite{lee_reaction_2007} derived analytical models, know as equimomental ellipsoids, to parameterize centroidal inertia. However, they introduce highly nonlinear constraints into the optimization, and discovering physically meaningful inertia trajectories solely through the ellipsoid parameterization is difficult.

\begin{figure}[tb!]
    \centering

    \includegraphics[width=\linewidth]{parkour_wide.pdf}
    \caption{Dynamic maneuvers that require high-fidelity control of angular momentum, such as ``parkour" (top) and cantering (bottom).}
    \label{fig:sim_ex}
    \vspace{-0.2in}
\end{figure}

Directly incorporating full-order models into a single-optimization is another common approach to complex whole-body motion generation. One option is to use full dynamics with either a hybrid \cite{hereid_3d_2016,mastalli_crocoddyl_2020,li_hybrid_2020} or contact-implicit \cite{posa_direct_2014,neunert_trajectory_2017,tassa_synthesis_2012,mordatch_discovery_2012} formulation. Alternatively, centroidal dynamics are utilized in \cite{dai_whole_2014,farshidian_efficient_2017,sleiman_constraint_2021} with full kinematics through centroidal momentum matrix (CMM) \cite{orin_centroidal_2013}. Notably, the differential dynamic programming (DDP)-based indirect approaches, such as FDDP \cite{mastalli_crocoddyl_2020} and constrained SLQ \cite{sleiman_constraint_2021}, have shown impressive results in terms of computational efficiency and constraint handling. Additional efforts have also been made in \cite{tassa_control_2014,xie_differential_2017,howell_altro_2019} to enforce general state, input equality and inequality constraints. However, a both efficient and general-purpose solver, capable of handling nonlinear full-order models, does not yet exist especially for highly agile acrobatic motions.

\begin{figure*}[t]
    \centering
    \smallskip
    \smallskip
    \smallskip
    \includegraphics[width=7in]{system_diagram_v3.pdf}
    \caption{Overall planning and control framework that consists of offline trajectory generation (red) and online execution (blue). The input to offline block is specifically decided by the desired motion. Full-body, centroidal, and contact information are passed from offline generation to online control.}
    \label{fig:system_diagram}
    \vspace{-0.2in}
\end{figure*}

By decoupling single-optimization into simpler sub-problems, it is still possible to leverage the benefits of both reduced and full order models. One strategy is through a hierarchical optimization. In \cite{carpentier_multicontact_2018,kudruss_optimal_2015}, a centroidal optimization solves for momentum and contact force trajectories, then inverse kinematics is applied to generate joint motions. Similarly, Nguyen et al. \cite{nguyen_contact_2021} proposed a pipeline, first solving an SRBM-based centroidal optimization, then using its solution as tracking objective for a full dynamics based whole-body optimization. However, additional constraints projecting whole-body information to the centroidal level are required to ensure feasible solutions. As pointed out by \cite{budhiraja_dynamics_2019}, these constraints are often computationally expensive and difficult to define, especially when involving angular momentum (AM). Another strategy would be to iteratively solve the two optimizations in an alternating fashion until the convergence of a common set of dynamics variables is achieved (dynamic consensus) \cite{herzog_structured_2016,ponton_efficient_2021,budhiraja_dynamics_2019,zhou_accelerated_2020}. However, existing frameworks compromise in either of the two optimizations: only full kinematics models (instead of dynamics ones) are considered in \cite{herzog_structured_2016,ponton_efficient_2021} during its whole-body optimization, while the centroidal optimizations proposed in \cite{budhiraja_dynamics_2019, zhou_accelerated_2020} did not explicitly optimize over foot placements. To our knowledge, no alternating optimization frameworks currently exist that both optimize over contact and use inertia-aware reduced order models at the centroidal level, while formulating constrained full dynamics DDP at the whole-body level.

In this work, we propose an alternating centroidal and whole-body optimization framework. Compared to existing alternating frameworks, we improve the centroidal optimization by incorporating the equimomental ellipsoid parameterization of centroidal inertia \cite{lee_reaction_2007} with SRBM, and leveraging the inertia feedback from whole-body optimization to enable a wide range of motions previously unachievable with standard centroidal optimization, such as ``inertia shaping" \cite{lee_reaction_2007} that aims to reach a goal pose through precise control of centroidal inertia. In addition, our whole-body optimization is designed to only incorporate minimal constraints to track nominal trajectories from the centroidal level in its objective; therefore, efficient methods such as DDP can be used. We ensure solution feasibility, in terms of constraint satisfaction, by handling most task-space constraints at the centroidal level via direct collocation, and improve solution quality at the whole-body level through dynamics consensus.

To validate the solution trajectory generated by our alternating optimization framework on real robots, we have developed a novel discrete-time finite-horizon model predictive tracking controller (MPC). Many predictive tracking controllers exist in the legged locomotion control literature. However, to achieve reasonable control frequency, it is common to use reduced order models within the problem formulation \cite{di_carlo_dynamic_2018, bledt_regularized_2020}. Although efficient and capable, these controllers are primarily designed for stabilizing fast periodic gaits, and often require expertly designed cost heuristics to generate desired motion \cite{bledt_regularized_2020}. These controllers also do not explicitly solve for joint variables which makes precisely tracking angular-momentum-rich joint motions impossible without introducing additional layers of whole-body kinematics/dynamics optimization \cite{kim_highly_2019, wang_unified_2021}. Our implementation of MPC aims to address these short comings. Specifically, the SRBM is augmented with centroidal momentum terms, and the joint variables are introduced through the CMM. With some mild assumptions, the proposed model can be linearized, which not only allows for a fast and efficient QP formulation of the MPC problem similar to \cite{di_carlo_dynamic_2018}, but also grants us direct control over the joint variables for direct momentum tracking.

Our core contributions are highlighted as follows:
\begin{itemize}
    \item We propose an alternating centroidal, whole-body optimization scheme targeting a wide range of motions. Inertia shaping tasks can be performed due to our inertia-aware SRBM based centroidal optimization. We also demonstrate improvement in solution feasibility compared to the standard whole-body optimization, via dynamics consensus.

    \item We formulate a convex model predictive controller through a novel linear transformation of the standard full centroidal dynamics.
    Our MPC is capable of tracking momentum-rich motions by jointly optimizing over contact forces and joint accelerations, while still retaining model simplicity for fast online computation.

    \item We show that the proposed framework generates high-fidelity contact and momentum-rich agile maneuvers, and demonstrate the capability of our controller on the MIT Mini-Cheetah platform.

\end{itemize}

\section{System Overview}\label{sec:system_overview}

We propose a hierarchical planning and control framework in this paper, as shown in Fig. \ref{fig:system_diagram}. The inputs to our framework consist of the desired motions, specified by the user (e.g. trot, canter, twisting jump, etc.), along with the associated contact sequence (not including contact locations), total motion duration, initial and final poses, and the terrain map. The offline trajectory planner takes these inputs, and solves centroidal and whole-body optimizations in an alternating fashion, until dynamics consensus is achieved. The resulting centroidal momentum, base motion, and joint trajectories are modified and tracked online using the proposed full centroidal MPC, which re-computes foot contact forces and joint accelerations based on state feedback from the robot. First of all, let us establish the system models.

\subsection{System Modeling}
Consider the standard floating base model of a legged robot, with an unactuated $6$-DoF base and a set of $n$-DoF fully-actuated limbs. The equations of motion is given by:
\begin{equation} \label{eq:fbEoM}
    \resizebox{0.65\hsize}{!}{$
    \begin{aligned}
        &\mathbf{M}(\mathbf{q}) \mathbf{\dot{v}}
        + \mathbf{C}(\mathbf{q, v})
        = \underbrace{ \begin{bmatrix} \textbf{0} \\ \textbf{I} \end{bmatrix} }_{\mathbf{B}} \vg \tau
        + \mathbf{J_{\rm c}^{\top}} \mathbf{F_{\rm c}} &&

    \end{aligned}
    $}

\end{equation}
where $\mathbf{q} = \begin{bmatrix} \mathbf{q_{\rm b}}^{\top} & \mathbf{q_{\rm j}}^{\top} \end{bmatrix}^{\top} \in \mathbb{R}^{n_{q}} \vspace{1mm}$, $\mathbf{v} = \begin{bmatrix} \vg{\nu}_{\rm b}^{\top} & \mathbf{\dot{q}_{\rm j}}^{\top} \end{bmatrix}^{\top} \in \mathbb{R}^{n_{v}}$ are the generalized coordinates and velocities partitioned in base and joint variables. The subscripts $\rm b$ and $\rm j$ are base and joint related quantities, respectively. Base coordinates $\mathbf{q_{\rm b}} = \begin{bmatrix} \mathbf{p_{\rm b}}^{\top} & \vg{\theta_{\rm b}}^{\top} \end{bmatrix}^{\top} \in \mathbb{R}^{6}$ are partitioned as base position and orientation. Base velocities $\vg{\nu}_{\rm b} = \begin{bmatrix} \mathbf{\dot{p}_{\rm b}}^{\top} & \vg{\omega}^{\top} \end{bmatrix}^{\top} \in \mathbb{R}^{6} \vspace{1mm}$ are partitioned as base linear and angular velocities in the world frame. $\mathbf{M(q)}$ is the joint space mass matrix. $\mathbf{C(q, v)}$ captures the nonlinear effects. $\vg{\tau} \in \mathbb{R}^{n_{j}}$ denotes joint torques. $\mathbf{J_{\rm c}(q)}$ is the stacked contact Jacobian and $\mathbf{F_{\rm c}}$ is the stacked contact reaction force vector. Without loss of generality, an implicit assumption $\mathbf{v} = \mathbf{\dot{q}}$ is made for the subsequent development of centroidal and whole-body optimization. This model will be referred to as the \textbf{full dynamics} of the robot.

Assuming full control authority over the joint variables $\mathbf{q_{\rm j}}$, and applying appropriate transformations, Eq. (\ref{eq:fbEoM}) can be converted into Newton-Euler equations of motion about the robot's center-of-mass (CoM) \cite{spong2006robot}:
\begin{equation}\label{eq:cenDyn}
  \resizebox{0.5\hsize}{!}{$
    \mathbf{\dot{h}} = \begin{bmatrix} \mathbf{\dot{k}} \\ \mathbf{\dot{l}} \end{bmatrix}
    = \begin{bmatrix} \sum\limits_{\rm j = 0}^{n_{f}-1}\mathbf{f_{\rm j}} + m \mathbf{g} \\ \sum\limits_{\rm j=0}^{n_{f}-1}(\mathbf{p_{\rm j}} - \mathbf{r}) \times \mathbf{f_{\rm j}} \end{bmatrix}
    $}
\end{equation}
where it is assumed that the robot has $n_f$ feet with point contacts. $\mathbf{h} = \begin{bmatrix} \mathbf{k}^{\top} & \mathbf{l}^{\top}\end{bmatrix}^{\top}$ are the linear and angular centroidal momentum. $\mathbf{p_{\rm{j}}}$ and $\mathbf{r}$ are the individual contact and CoM position in the world frame. $\mathbf{f_{\rm j}}$ are the individual foot contact forces. This reduced-order model will be referred to as the \textbf{centroidal dynamics} throughout this paper. Unless otherwise stated, all vector quantities in this paper are expressed in the world frame.

\section{Offline Trajectory Generation}\label{sec:trajectoryGeneration}

We present an alternating centroidal, whole-body trajectory optimization scheme in this section. Detailed formulations for individual components are provided in Sec. \ref{subsec:cen_traj_gen} and \ref{subsec:wbd_traj_gen}. The alternating update scheme and the importance of dynamics consensus are discussed in Sec. \ref{subsec:alternating_updates}.

\subsection{Centroidal Trajectory Optimization}\label{subsec:cen_traj_gen}
The TO problem to be solved at the centroidal level utilizes the \textbf{centroidal dynamics} in Eq. (\ref{eq:cenDyn}). We transcribe the optimization problem into a nonlinear program through direct collocation, due to its natural capability of constraint handling.
Total duration $T$ of the desired motion is equally divided into $N$ knot points. Contact sequence and timing per foot are also predefined. Formulation \ref{formulation:OCP for centroidal} details the centroidal optimization.

The centroidal decision set $\vg \phi_{\rm cen}$ consists of CoM position $\v r$, CoM linear velocity $\dot{\v r}$, angular excursion \cite{zordan_control_2014} with Euler angle parameterization $\vg \theta$ and its rate $\dot{\vg \theta}$, centroidal angular momentum $\v l$ and its rate $\dot{\v l}$, $j^{\rm th}$ foot's position $\v p_{\rm j}$ and contact force $\v f_{\rm j}$.
$\v e$ and $\vg \gamma$ are inertia variables characterizing the principle semi-axes and orientation of the equimomental ellipsoid \cite{lee_reaction_2007} (see the ellipsoid illustration in Fig.~\ref{fig:TO_system_diagram}). Ellipsoid orientation is again parameterized using Euler angles. All components of $\vg \phi_{\rm cen}$ belongs to $\mathbb{R}^3$.

The optimization objective contains a user-defined cost $\mathcal{L}_{\rm cen}(\cdot)$ which specifies various heuristics, such as penalizing the velocity of foot position to discourage violent movements. In addition, a tracking cost $\Psi_{\rm cen}(\vg \phi_{\rm cen},\vg \phi_{\rm wbd}^{\rm ref})$ that minimizes the tracking error to the whole-body reference trajectories (solution to Formulation \ref{formulation:OCP for whole-body} and the variable $\vg \phi_{\rm wbd}^{\rm ref}$ is defined in Formulation \ref{formulation:OCP for centroidal}) are included. Benefits of this tracking cost will be discussed in Sec.~\ref{subsec:alternating_updates}.
\begin{equation*}
  \resizebox{0.75\hsize}{!}{$
\begin{aligned}
     \Psi_{\rm cen} = &\|\v r - \v r^{\rm ref}_{\rm wbd} \|^2_{\mathbf{Q^{\rm cen}_{\rm r}}}+\|m \mathbf{\dot{r}} - \mathbf{k^{\rm ref}_{\rm wbd}} \|^2_{\mathbf{Q^{\rm cen}_{\rm k}}}+\\
     &\|\v l - \v l^{\rm ref}_{\rm wbd} \|^2_{\mathbf{Q^{\rm cen}_{\rm l}}}+\|\v I_{\rm ellip} - \v I_{\rm crb} \|^2_{\mathbf{Q^{\rm cen}_{\rm I}}}+\\&\sum_{j=0}^{n_f-1}(\|\v p_j - \v p_{j,\rm wbd}^{\rm ref} \|^2_{\mathbf{Q^{\rm cen}_{\rm p}}})
\end{aligned}$}
\end{equation*}
Symmetric positive definite weighting matrices for tracking CoM position, linear momentum (LM), angular momentum (AM), and foot positions are denoted as $\mathbf{Q^{\rm cen}_{\rm r}}$, $\mathbf{Q^{\rm cen}_{\rm k}}$, $\mathbf{Q^{\rm cen}_{\rm l}}$, and $\mathbf{Q^{\rm cen}_{\rm p}}$, respectively. Inertia tracking is achieved through minimizing the $\mathbf{Q^{\rm cen}_{\rm I}}$-weighted Frobenius norm of the difference between the inertia tensors of the equimomental ellipsoid $\mathbf{I_{\rm ellip}}$ and that of the reference composite rigid body (CRB) $\mathbf{I_{\rm crb}}$, both expressed in the world frame. Centroidal optimization includes the following constraints:

$\mathbf{Dynamic \ constraint}$:
Without any simplification, the discrete-time non-linear \textbf{centroidal dynamics} are used. Eq.~(\ref{eq:dynamics}) explicitly includes the relationship between centroidal angular momentum, ellipsoid inertia tensor, and CoM angular velocity. The inertia tensor $\v I_{\rm ellip}$ in the world frame is constructed from the rotational transformation of the inertia tensor in the base frame, $\mathbf{I_{\rm ellip}} = \mathbf{R}(\vg \gamma) _{\mathcal{B}}\v I_{\rm ellip}(\mathbf{e}) \mathbf{R}(\vg \gamma)^{\top}$ with $\mathbf{R}(\vg \gamma)$ being the rotation transformation between the ellipsoid base frame and world frame. Inertia tensor in the base frame is defined as $_{\mathcal{B}}\v I_{\rm ellip} = diag([\mathbf{I}^{\rm xx}, \mathbf{I}^{\rm yy}, \mathbf{I}^{\rm zz}]^T)$  with $\mathbf{I}^{xx} = \frac{1}{5}m(\v {e^{\rm y}}^2 + \v {e^{\rm z}}^2)$, $\mathbf{I}^{yy} = \frac{1}{5}m(\v {e^{\rm x}}^2 + \v {e^{\rm z}}^2)$, and $\mathbf{I}^{zz} = \frac{1}{5}m(\v {e^{\rm x}}^2 + \v {e^{\rm y}}^2)$.

$\mathbf{Ellipsoid \ mass \ constraint}$: To ensure the equimomental ellipsoid has the same mass as the robot, ellipsoid semi-axes, $\mathbf{e}$, must be constrained. Ellipsoid mass $m = \frac{4}{3} \pi \mathbf{e^{\rm x}} \mathbf{e^{\rm y}} \mathbf{e^{\rm z}} \rho$ with $\rho$ being the average mass density. $\rho$ can be obtained from the reference CRB inertia tensor $\mathbf{I_{\rm crb}}$ for each time step \cite{orin_centroidal_2013}. When the CRB inertia reference is unavailable, $\rho$ is calculated using the initial robot configuration and kept constant.

$\mathbf{Integration \ constraint}$:
By adopting trapezoidal integration scheme, the quantities $\v r$, $\dot{\v r}$, $\vg \theta$ and $\v l$ are approximated as piecewise linear functions as shown in Eqs.~(\ref{eq:integration_r}) - (\ref{eq:integration_L}). The final knot points are set to be free but subject to final conditions if needed. $\Delta t$ denotes one time step.

$\mathbf{Frictional \ constraint}$:
Friction cone constraint and unilateral constraint on the contact forces are both modeled in Eqs.~(\ref{eq:force_non_zero}) - (\ref{eq:force_zero}). Specifically, friction cones are approximated as friction pyramids $\mathcal{F}$ to remove nonlinearity. The contact normal forces are non-negative, and always equal to zero during the swing phase.

$\mathbf{Terrain \ constraint}$:
For each foot during the stance phase, Eq.~(\ref{eq:terrain_1}) defines the no-slip condition. While Eq.~(\ref{eq:terrain_2}) guarantees that foot placement remains on the terrain surface, as specified by the height map $h_{\rm terrain}$ (see \cite{winkler_gait_2018} for details).

$\mathbf{Kinematic \ constraint}$:
Foot position range-of-motion constraints, Eq.~(\ref{eq:rom}), are defined, to heuristically bound each foot's movement within a fixed sized box $\mathcal{R}$.

\begin{subfigures}
\begin{figure}
    \vspace{-0.5mm}
\end{figure}
\begin{formulation}[t]
\caption{Centroidal trajectory optimization}\label{formulation:OCP for centroidal}
\begin{subequations}
\resizebox{1.0\linewidth}{!}{
  \begin{minipage}{\linewidth}
    \begin{align}
        \nonumber &\hspace{2.1cm} \underset{\vg \phi_{\rm cen}}{\min} \ \mathcal{L}_{\rm cen}(\vg \phi_{\rm cen}) + \Psi_{\rm cen}(\vg \phi_{\rm cen},\vg \phi_{\rm wbd}^{\rm ref})
        \\\nonumber
        &(\textit{\rm Variables}) \hspace{0.5cm} \vg \phi_{\rm cen}[i] = [\v r[i]^{\top}, \dot{\v r}[i]^{\top}, \vg \theta[i]^{\top}, \dot{\vg \theta}[i]^{\top},\v l[i]^{\top}, \\\nonumber
        &\hspace{3.5cm}  \dot{\v l}[i]^{\top}, \v p_j[i]^{\top}, \v f_j[i]^{\top},\underbrace{\v e[i]^{\top}, \vg \gamma[i]^{\top}}_{\text{inertia vectors}}
        ]^{\top}\\\nonumber & \hspace{2.1cm} \vg \phi_{\rm wbd}^{\rm ref}[i]=[\v r^{\rm ref}_{\rm wbd}[i]^{\top},\mathbf{\dot{r}^{\rm ref}_{\rm wbd}}[i]^{\top},\v l^{\rm ref}_{\rm wbd}[i]^{\top},\\
        \nonumber & \hspace{3.7cm} \v p^{\rm ref}_{j, \rm wbd}[i]^{\top}, \v I_{\rm crb}[i]^{\top}]^{\top}\\\nonumber
        &\hspace{2.1cm} \forall i = 0, \ldots, N-1, j=0,\ldots,n_f-1\\\label{eq:dynamics}
        &(\textit{\rm Dynamics}) \ \text{s.t.}  \begin{bmatrix} m \ddot{\v r}[i]\\\
        \dot{\v l}[i]\\\
        \v l[i]
        \end{bmatrix}= \begin{bmatrix}
        \sum\nolimits_{j} \v f_j[i] + m \v g \\
	    \sum\nolimits_j (\v{p}_j[i] - \v r[i]) \times \v f_j[i] \\
	    \v I_{\rm ellip}(\v e[i], \vg \gamma[i])\vg \omega[i]
	    \end{bmatrix} \\\label{eq:ellipsoidMass}
	    &(\textit{\rm Ellipsoid}) \hspace{0.52cm} m = \frac{4}{3} \pi \mathbf{e^{\rm x}}[i] \mathbf{e^{\rm y}}[i] \mathbf{e^{\rm z}}[i] \rho[i] \\\label{eq:integration_r}
        &(\textit{\rm Integration}) \hspace{0.3cm} \v r[i+1] - \v r[i] = \frac{\Delta t}{2}(\dot{\v r}[i+1]+\dot{\v r}[i])\\\label{eq:integration_r_dot}
        &\hspace{2.1cm} \dot{\v r}[i+1] - \dot{\v r}[i] = \frac{\Delta t}{2}(\ddot{\v r}[i+1]+\ddot{\v r}[i])\\\label{eq:integration_theta}
        &\hspace{2.1cm} \vg \theta[i+1] - \vg \theta[i] = \frac{\Delta t}{2}(\dot{\vg \theta}[i+1]+\dot{\vg \theta}[i])\\\label{eq:integration_L}
        &\hspace{2.1cm} \v l[i+1] - \v l[i] = \frac{\Delta t}{2}(\dot{\v l}[i+1]+\dot{\v l}[i])\\\label{eq:force_non_zero}
        &(\textit{\rm Frictional})

    \hspace{0.45cm} \v f_j[i] \cdot \v n(\v p_j^{xy}[i]) \ge 0, \ \forall i \in \mathcal{C}_j\\\label{eq:force_friction_cone}
    &\hspace{2.1cm} \v f_j[i] \in \mathcal{F}(\mu, \v n, \v p_j^{xy}[i]), \ \forall i \in \mathcal{C}_j\\\label{eq:force_zero}

     &\hspace{2.1cm} \v f_j[i] = 0, \ \forall i \notin \mathcal{C}_j\\\label{eq:terrain_1}
        &(\textit{\rm Terrain})

    \hspace{0.82cm} \v p_j^{xy}[i+1] - \v p_j^{xy}[i]=0, \ \forall i \in \mathcal{C}_j\\\label{eq:terrain_2}
    &\hspace{2.12cm} \v p_j^z[i]=h_{\rm terrain}(\v p_j^{xy}[i]), \ \forall i \in \mathcal{C}_j\\\label{eq:rom}
        &(\textit{\rm Kinematics}) \hspace{0.22cm} \v p_j[i] \in \mathcal{R}_j(\v r[i],\vg \theta[i])
    \end{align}
     \end{minipage}
     }
\end{subequations}
\end{formulation}
\end{subfigures}

\subsection{Whole-Body Trajectory Optimization}\label{subsec:wbd_traj_gen}
For legged systems, leveraging full-body dynamics, such as the angular momentum of limbs, is critical for generating dynamic motions. Since the contact sequence is provided beforehand, the whole-body optimization is formulated in a hybrid dynamics fashion. The objective is to track quantities computed from centroidal optimization. Formulation \ref{formulation:OCP for whole-body} details this whole-body optimization.

The decision set $\vg \phi_{\rm wbd}$ consists of the generalized position $\v q$ and velocity $\dot{\v q}$, joint torques $\vg \tau$ and contact forces $\mathbf{F_{\rm c}}$. Similar to the centroidal optimization, the objective includes a user-defined heuristic cost $\mathcal{L}_{\rm wbd}$, such as torque minimization, and a tracking cost $\Psi_{\rm wbd}(\vg \phi_{\rm wbd},\vg \phi_{\rm cen}^{\rm ref})$ to minimize the deviation from the centroidal reference trajectories (solution to Formulation \ref{formulation:OCP for centroidal}).
\begin{equation*}
  \resizebox{0.85\hsize}{!}{$
    \begin{aligned}
      \Psi_{\rm wbd}=& \|\mathcal{X}_r(\v q) - \v r^{\rm ref}_{\rm cen} \|^2_{\mathbf{Q^{\rm wbd}_{\rm r}}}+\|\v A(\v q)\dot{\v q} - \v h^{\rm ref}_{\rm cen}\|^2_{\mathbf{Q^{\rm wbd}_{\rm h}}}+\\&\sum_{j=0}^{n_f}(\|\mathcal{X}_j(\v q) - \v p^{\rm ref}_{j,\rm cen}\|^2_{\mathbf{Q^{\rm wbd}_{\rm p}}})
    \end{aligned}$}
\end{equation*}
where the quadratic weighting matrices for tracking CoM position/orientation, centroidal momentum, and foot positions are denoted as $\mathbf{Q^{\rm wbd}_{\rm r}}$, $\mathbf{Q^{\rm wbd}_{\rm h}}$, and $\mathbf{Q^{\rm wbd}_{\rm p}}$, respectively. $\mathcal{X}_r(\cdot)$ and $\mathcal{X}_j(\cdot)$ map the generalized position to CoM position and the $j \rm th$ end-effector (EE)'s position. The centroidal momentum matrix, $\mathbf{A(q)} \in \mathbb{R}^{6 \times (6 + n_j)}$, referring to Eq. (\ref{eq:cmmToGeneralizedVel}), provides a linear relationship between the centoridal momentum $\mathbf{h}$ and generalized velocities $\mathbf{\dot{q}}$ \cite{orin_centroidal_2013}. The hybrid dynamic constraint includes both the non-impact \textbf{full dynamics} Eq. (\ref{eq:fbEoM}), and the impact phases described in Eq.~(\ref{eq:hybrid_dynamics}), where $\vg \Lambda$ stands for the contact impulse, and $\dot{\v q}^-$ and $\dot{\v q}^+$ are instantaneous joint velocities before and after the impact. Additionally, a no-slip contact constraint (\ref{eq:no_slip}) is added for the stance foot. Joint limit, torque limit, and friction cone are also modeled as shown in Eqs. (\ref{eq:joint_torque_limits}) - (\ref{eq:wbd_friction_cone}).

\begin{subfigures}
\begin{figure}
    \vspace{0mm}
\end{figure}
\begin{formulation}[t]
\caption{Whole-body trajectory optimization}
\label{formulation:OCP for whole-body}
\begin{subequations}
\resizebox{1.0\linewidth}{!}{
  \begin{minipage}{\linewidth}
    \begin{align}
        \nonumber &\hspace{1.95cm} \underset{\vg \phi_{\rm wbd}}{\min} \ \mathcal{L}_{\rm wbd}(\vg \phi_{\rm wbd}) + \Psi_{\rm wbd}(\vg \phi_{\rm wbd},\vg \phi_{\rm cen}^{\rm ref})\\\nonumber
        &(\textit{\rm Variables}) \quad \vg \phi_{\rm wbd}[i] = [\v q[i]^{\top}, \dot{\v q}[i]^{\top}, \vg \tau[i]^{\top}, \mathbf{F_{\rm c}}[i]^{\top}]^{\top}\\\nonumber & \hspace{1.95cm} \vg \phi_{\rm cen}^{\rm ref}[i]=[\v r^{\rm ref}_{\rm cen}[i]^{\top},\v h^{\rm ref}_{\rm cen}[i]^{\top},\v p^{\rm ref}_{j,\rm cen}[i]^{\top}]^{\top}\\\nonumber
        &\hspace{1.95cm} \forall i = 0, 1, \ldots, N-1, j=0,\ldots,n_f-1\\\label{eq:hybrid_dynamics}
        &(\textit{\rm Dynamics}) \quad         \text{s.t.} \begin{cases}
        \v M(\v q) \v {\ddot q} + \v C(\v q, \v {\dot q})  = \v B \vg \tau + \mathbf{J_{\rm c}^{\top}} \mathbf{F_{\rm c}}\\
        \v M(\v q)\dot{\v q}^+ - \v M(\v q)\dot{\v q}^- = \mathbf{J_{\rm c}^{\top}} \vg \Lambda
        \end{cases}\\\label{eq:no_slip}
        &(\textit{\rm No-slip}) \hspace{1.45cm} \v J_c \ddot{\v q}+\mathbf{\dot{J}_{\rm c}}\dot{\v q} = 0\\\label{eq:joint_torque_limits}
        &(\textit{\rm Limits}) \hspace{1.6cm} \v q \in \mathcal{J}, \vg \tau \in \mathcal{T}\\\label{eq:wbd_friction_cone}
        &(\textit{\rm Friction cone}) \hspace{0.6cm} \vg \lambda_c \in \mathcal{K}
    \end{align}
    \end{minipage}
}
\end{subequations}
\end{formulation}
\end{subfigures}

Considering the computational efficiency, we then choose the DDP-based method implemented in \cite{mastalli_crocoddyl_2020} to solve the whole-body optimization. Briefly, the dynamic and no-slip constraints are coupled to derive the joint acceleration and contact forces simultaneously, so that the forward pass for DDP can be performed. The joint limit, torque limit and friction cone constraints are encoded as quadratic barrier functions inside the objective function. Readers are referred to \cite{mastalli_crocoddyl_2020} for more details. Based on the above whole-body trajectory optimization, we are able to generate versatile whole-body motions through tracking reference trajectories provided by the centroidal optimization.

\subsection{Alternating CEN-WBD Optimization}\label{subsec:alternating_updates}
We perform an alternating update scheme based on the aforementioned centroidal optimization (CEN) and whole-body optimization (WBD) problems. Fig.~\ref{fig:TO_system_diagram} demonstrates the proposed updating sequence, where the CEN optimization is first solved and then followed by the WBD optimization. The reference trajectories are updated accordingly after each CEN or WBD optimization solve. $\mathcal{I}_{\rm crb}$ maps the generalized coordinates to the CRB inertia tensor $\v I_{\rm crb}$. Similar to \cite{herzog_structured_2016,ponton_efficient_2021,budhiraja_dynamics_2019,zhou_accelerated_2020}, this process, referred to as one alternating iteration, is iterated until the dynamics consensus between CEN and WBD is achieved. Moreover, our proposed consensus quantities include not only CoM position, centroidal momentum, and EE locations, but also the equimomental inertia tensors ($\v I_{\rm ellip}$ in CEN and $\v I_{\rm crb}$ in WBD).  This is achieved by enabling the variables, constraints, and costs that involve ellipsoid inertia vectors $\v e$ and $\v v$ in CEN, which allows generating robot behaviors that cannot be captured by only using centroidal dynamics such as zero-gravity reorientation (see Remark \ref{remark_1}). For the other motions unnecessary to parameterize centroidal inertia, we disable the ellipsoid related terms in CEN and directly use $\v I_{\rm crb}$ from WBD. This hybrid formulation reduces the complexity of directly solving Formulation \ref{formulation:OCP for centroidal} and saves the solve time. More results will be shown in Sec.~\ref{sec:experiments}.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=3.4in]{TO_diagram_v4.pdf}
    \caption{The block diagram that demonstrates the proposed alternating centroidal, whole-body optimization. }
    \label{fig:TO_system_diagram}
    \vspace{-0.2in}
\end{figure}

\begin{remark}
\label{remark_1}
We introduce two sets of orientation representations in CEN: angular excursion $\vg \theta$ and equimomental ellipsoid orientation $\v v$. The former one can be seen as the angular analog of CoM and only changed under external torques \cite{zordan_control_2014}, while the latter one can be affected given any limb movements even without external forces \cite{lee_reaction_2007}.
\end{remark}

\section{Full Centroidal Convex MPC}
This section describes the mathematical formulation of our proposed tracking controller. Starting with the derivation of the linearized full centroidal dynamics, we present a QP-based convex-MPC problem. We then discuss the benefits of our controller and provides the comparison with the original convex-MPC \cite{di_carlo_dynamic_2018} in Sec. \ref{sec:comparisonOriginalMPC}.

\subsection{Full Centroidal Dynamics}

Standard full centroidal dynamics as shown in Eq. (\ref{eq:cenDyn}) has been widely used in the legged systems literature \cite{dai_whole_2014, farshidian_efficient_2017}. We propose a linearized variant of this model, and demonstrate its benefits within the context of MPC. To make this linear model unambiguous, the transformation between the time derivative of generalized base coordinates and the generalized base velocities is explicitly stated here:
\begin{equation}\label{eq:quatRateAnglVel}
  \resizebox{0.4\hsize}{!}{$
    \mathbf{\dot{q}_{\rm b}} = \underbrace{\begin{bmatrix} \textbf{I} & \textbf{0}_{3 \times 3} \\ \textbf{0}_{4 \times 3} & \mathbf{W}(\vg{\theta}_{\rm b}) \end{bmatrix}}_{\mathbf{T_{\rm b}}(\mathbf{q})} \vg{\nu}_{\rm b}
    $}
\end{equation}
where $\mathbf{W(\vg \theta_{\rm b})}$ represents the velocity transform for different choices of orientation parameterization, see \cite{robotCheatSheet}.

Using the CMM, $\mathbf{A(q)} \in \mathbb{R}^{6 \times (6 + n_j)}$, to relate centroidal momentum to generalized velocities:
\begin{equation}\label{eq:cmmToGeneralizedVel}
  \resizebox{0.4\hsize}{!}{$
    \mathbf{h} = \underbrace{\begin{bmatrix} \mathbf{A_{\rm b}}(\mathbf{q}) & \mathbf{A_{\rm j}}(\mathbf{q}) \end{bmatrix}}_{\mathbf{A}(\mathbf{q})} \begin{bmatrix} \vg{\nu}_{\rm b} \\ \mathbf{\dot{q}_{\rm j}} \end{bmatrix}
$}
\end{equation}
Rearranging Eq. (\ref{eq:cmmToGeneralizedVel}) for $\vg{\nu_{\rm b}}$ and using Eq. (\ref{eq:quatRateAnglVel}), one has
\begin{equation}\label{eq:baseCoordRate}
    \mathbf{\dot{q}_{\rm b}} = [~\underbrace{\mathbf{T_{\rm b}}(\mathbf{q}) \mathbf{A_{\rm b}^{-1}}(\mathbf{q})}_{\mathbf{\widetilde{A_{\rm h}}}(\mathbf{q})}~~~\underbrace{-\mathbf{T_{\rm b}}(\mathbf{q}) \mathbf{A_{\rm b}^{-1}}(\mathbf{q}) \mathbf{A_{\rm j}}(\mathbf{q})}_{\mathbf{\widetilde{A_{\rm j}}}(\mathbf{q})}~] \begin{bmatrix} \mathbf{h} \\ \mathbf{\dot{q}_{\rm j}} \end{bmatrix}
\end{equation}
where $\mathbf{\widetilde{A_{\rm h}}}(\mathbf{q})$ and $\mathbf{\widetilde{A_{\rm j}}}(\mathbf{q})$ will be referred to as transformed CMMs. Defining an augmented state $\mathbf{x} = \begin{bmatrix} \mathbf{h}^{\top},\mathbf{q_{\rm b}}^{\top},\mathbf{q_{\rm j}}^{\top},\mathbf{\dot{q}_{\rm j}}^{\top},\mathbf{g}^{\top} \end{bmatrix}^{\top} \vspace{1mm}$, and combining Eqs. (\ref{eq:cenDyn}) and (\ref{eq:baseCoordRate}), the augmented state dynamics can be written as:
\begin{equation}\label{eq:fullCenDyn}
  \resizebox{0.85\hsize}{!}{$
  \begin{aligned}
    \begin{bmatrix} \mathbf{\dot{h}} \\ \mathbf{\dot{q}_{\rm b}} \\ \mathbf{\dot{q}_{\rm j}} \\ \mathbf{\ddot{q}_{\rm j}} \\ \mathbf{\dot{g}} \end{bmatrix}
    &= \underbrace{\begin{bmatrix} \textbf{0} & \textbf{0} & \textbf{0} & \textbf{0} & m \begin{bmatrix} \textbf{I} \\ \textbf{0} \end{bmatrix} \\
    \mathbf{\widetilde{A_{\rm h}}}(\mathbf{q}) & \textbf{0} & \textbf{0} & \mathbf{\widetilde{A_{\rm j}}}(\mathbf{q}) & \textbf{0} \\ \textbf{0} & \textbf{0} & \textbf{0} & \textbf{I} & \textbf{0} \\ \textbf{0} & \textbf{0} & \textbf{0} & \textbf{0} & \textbf{0} \\ \textbf{0} & \textbf{0} & \textbf{0} & \textbf{0} & \textbf{0} \end{bmatrix}}_{\mathbf{H}(\mathbf{q})}
      \begin{bmatrix} \mathbf{h} \\ \mathbf{q_{\rm b}} \\ \mathbf{q_{\rm j}} \\ \mathbf{\dot{q}_{\rm j}} \\ \mathbf{g} \end{bmatrix} &&\\\
    &+  \underbrace{\begin{bmatrix} \textbf{I} & \cdots & \textbf{I} & \textbf{0} \\ \widehat{(\mathbf{p_{\rm 0}}-\mathbf{r})} & \cdots & \widehat{(\mathbf{p_{\rm n_f - 1}}-\mathbf{r})} & \textbf{0} \\ \textbf{0} & \cdots & \textbf{0} & \textbf{0} \\
    \textbf{0} & \cdots & \textbf{0} & \textbf{0} \\
    \textbf{0} & \cdots & \textbf{0} & \textbf{I} \\
    \textbf{0} & \cdots & \textbf{0} & \textbf{0}
    \end{bmatrix}}_{\mathbf{G}(\mathbf{p}, \mathbf{r})}
    \underbrace{
    \begin{bmatrix}
        \mathbf{f_{\rm 0}} \\
        \vdots \\
        \mathbf{f_{\rm n_f - 1}} \\
        \mathbf{a}
    \end{bmatrix}
    }_{\mathbf{u}}
  \end{aligned}
  $}
\end{equation}
We refer to this system as the \textbf{full centroidal dynamics}. The $\widehat{\cdot}$ operator transforms the cross-product operation into a matrix multiplication. Note that, $\mathbf{H}$ only depends on the joint configuration $\mathbf{q}$, and $\mathbf{G}$ only depends on contact and CoM locations, $\mathbf{p}$ and $\mathbf{r}$. If the robot follows the desired reference trajectory generated from our offline optimization, then foot placements, CoM location, and transformed CMMs from offline trajectory can be directly substituted into Eq. (\ref{eq:fullCenDyn}).
\begin{equation} \label{eq:LTVfullCenDyn}
  \resizebox{0.5\hsize}{!}{$
    \mathbf{\dot{x}} = \mathbf{H(q^{\rm ref})} \mathbf{x} + \mathbf{G(p^{\rm ref}, r^{\rm ref})} \mathbf{u}
$}
\end{equation}
which simplifies to a linear time-varying (LTV) system in the augmented state $\mathbf{x}$, with contact forces and joint accelerations as control inputs. This linearized dynamics will be used in the subsequent QP-based MPC formulation.

\subsection{Full Centroidal Convex MPC} \label{sec:fullCenConvexMPC}

Continuous LTV dynamics Eq. (\ref{eq:LTVfullCenDyn}) can be discretized. The resulting LTV discrete dynamics are given as:
\begin{equation} \label{eq:DiscreteLTVfullCenDyn}
    \mathbf{x}\scalebox{0.8}{$[i+1]$} = \mathbf{\widetilde{H}}\scalebox{0.8}{$[i]$} \mathbf{x}\scalebox{0.8}{$[i]$} + \mathbf{\widetilde{G}}\scalebox{0.8}{$[i]$} \mathbf{u}\scalebox{0.8}{$[i]$}
\end{equation}
where $\mathbf{\widetilde{H}}$ and $\mathbf{\widetilde{G}}$ are matrix exponential of $\mathbf{H}$ and $\mathbf{G}$ from Eq. (\ref{eq:LTVfullCenDyn}). Similar to the original convex-MPC \cite{di_carlo_dynamic_2018}, to reduce the size of the proposed MPC problem, state variables can be eliminated from the decision set by adopting a condensed dynamics formulation \cite{jerez_denseQP_2011}. At time step $k$, assuming a MPC prediction horizon of $n \in \mathbb{N}^{+}$ and aggregating Eq. (\ref{eq:DiscreteLTVfullCenDyn}) from $i = k$ to $i = k+n-1$ produces:
\begin{equation}
    \Vec{\mathbf{x}} = \mathbf{\widetilde{H^{\rm \mathcal{A}}}} \mathbf{x}[k] + \mathbf{\widetilde{G^{\rm \mathcal{A}}}} \Vec{\mathbf{u}}
\end{equation}
with aggregated states, \scalebox{0.8}{}$\Vec{\mathbf{x}} = \scalebox{1.3}{[} \mathbf{x}[k+1]^{\top} ~~ \mathbf{x}[k+2]^{\top} ~ \cdots ~ \mathbf{x}[k+n]^{\top} \scalebox{1.3}{]}^{\top}$. Weighting matrix $\mathbf{\widetilde{H^{\mathcal{A}}}}$ is given as
\begin{equation*}
    \resizebox{0.75\hsize}{!}{$
    \mathbf{\widetilde{H^{\mathcal{A}}}} = \begin{bmatrix} \prod\limits_{i = k}^{k}\mathbf{\widetilde{H}}[i]^{\top} & \prod\limits_{i = k}^{k+1}\mathbf{\widetilde{H}}[i]^{\top} & \cdots & \prod\limits_{i = k}^{k+n-1}\mathbf{\widetilde{H}}[i]^{\top} \end{bmatrix}^{\top}
    $}
\end{equation*}
and actuation matrix $\mathbf{\widetilde{G^{\mathcal{A}}}}$ is given as
\begin{equation*}
    \resizebox{0.9\hsize}{!}{$
    \mathbf{\widetilde{G^{\mathcal{A}}}} = \begin{bmatrix} \mathbf{\widetilde{G}}[k]\\ \prod\limits_{i = k+1}^{k+1}\mathbf{\widetilde{H}}[i] \mathbf{\widetilde{G}}[k] & \mathbf{\widetilde{G}}[k+1] & &  \text{\LARGE 0} \\ \vdots & \vdots & \ddots \\ \prod\limits_{i = k+1}^{k+n-1}\mathbf{\widetilde{H}}[i] \mathbf{\widetilde{G}}[k] & \prod\limits_{i = k+2}^{k+n-1}\mathbf{\widetilde{H}}[i] \mathbf{\widetilde{G}}[k+1] & \cdots & \mathbf{\widetilde{G}}[k+n-1] \end{bmatrix}
    $}
\end{equation*}
The decision variables for this dense QP formulation are aggregated control vectors $\Vec{\mathbf{u}}$ over the $n$-step prediction horizon, $\Vec{\mathbf{u}} = \scalebox{1.3}{[} \mathbf{u}[k]^{\top} ~~ \mathbf{u}[k+1]^{\top} ~ \cdots ~ \mathbf{u}[k+n-1]^{\top} \scalebox{1.3}{]}^{\top}$, where each control vector $\mathbf{u}$ consists of the foot contact forces $\mathbf{f}_{j}$ and the joint accelerations $\mathbf{a}$ for each actuated joint. It is straightforward to define the QP objective as a sum of tracking cost and a regularization term:
\begin{equation*}
    \vg \Psi_{\rm qp} = \| \mathbf{\widetilde{H^{\mathcal{A}}}}  \mathbf{x}[k] + \mathbf{\widetilde{G^{\mathcal{A}}}} \Vec{\mathbf{u}} - \Vec{\mathbf{x}}^{\rm{ref}} \|_{\mathbf{Q_{\rm{qp}}}} + \| \Vec{\mathbf{u}} \|_{\mathbf{R_{\rm{qp}}}}
\end{equation*}
where $\Vec{\mathbf{x}}^{\rm{ref}} = [\mathbf{h^{\rm ref^{\top}}}, \mathbf{q_{\rm b}^{\rm ref^{\top}}}, \mathbf{q_{\rm j}^{\rm ref^{\top}}}, \mathbf{\dot{q}^{\rm ref^{\top}}_{\rm j}}]^{\top}$ are given by the desired reference, generated through offline optimization, for the entire horizon $i = k+1$ to $i = k+n$. Friction pyramid constraint, joint acceleration constraint, and normal force constraint are also added into the QP formulation to ensure physically feasible contact forces and joint accelerations. The computed contact forces and joint accelerations from the proposed MPC-QP problem are fed into the inverse dynamics (ID) controller to generate joint torques as shown in Fig.~\ref{fig:system_diagram}.

\subsection{Comparison with the original Convex-MPC}\label{sec:comparisonOriginalMPC}
Compared to the original Convex-MPC \cite{di_carlo_dynamic_2018}, our formulation added the ability to directly track centroidal momentum quantities. This offers a tighter integration with our offline trajectory generation, and fully utilizes the generated reference linear and angular momentum trajectories to faithfully execute momentum-rich motions, such as twisting jump and zero-gravity body reorientation (see Sec. \ref{sec:experiments}). Unlike the original convex-MPC controller \cite{di_carlo_dynamic_2018}, where cost heuristics must be designed to generate such motions \cite{bledt_regularized_2020}, our controller allows momentum to be designed as part of the desired reference trajectory, and directly tracked online. In addition, leveraging CMMs in the system dynamics naturally allows the inclusion of joint variables as decision variables, thus eliminating the need for any additional kinematics optimizations commonly found in such control hierarchies \cite{kim_highly_2019, wang_unified_2021}. Another benefit when compared to the original is that our formulation does not assume zero roll direction rotation of the robot. On the downside, the ability to track momentum and joint trajectories directly increases the QP problem size, and our formulation is indeed slower to compute compared to the original, for the same length prediction horizon. However, by trading off between update frequency and prediction horizon, we still managed to achieve satisfactory tracking performance on real hardware (see Sec. \ref{sec:hardwareResults}).

The main assumption in our formulation is that the robot must follow the reference trajectory generated by offline planning for the approximation in Eq. (\ref{eq:LTVfullCenDyn}) to hold true. This assumption is reasonable because when the robot deviates from the offline reference, our MPC controller will generate new contact forces and joint accelerations, at a fast rate (see Sec. \ref{sec:hardwareResults}), to bring the robot state back to the reference.

\section{Results}\label{sec:experiments}
We perform various experiments both in simulation and hardware using the alternating CEN-WBD and MPC framework. Both the Unitree A1 and MIT Mini Cheetah \cite{katz_mini_2019} were used to verify the proposed methods on quadrupedal platforms. In simulation, dynamic maneuvers such as cantering, ``parkour", and zero gravity reorientation are shown to have stable and converging solutions. In addition, hardware tests were conducted with the Mini Cheetah to verify both the offline optimization and online MPC in the real world.
\subsection{Trajectory Generation and Simulation}

Example dynamic maneuvers generated using our offline alternating CEN-WBD optimization are presented here. We demonstrate the solution quality of these maneuvers either in visualization or physics simulation, and discuss the importance of dynamics consensus. Aggregated solve times and residuals for all the desired motions are documented in Table \ref{tab:solving_time_and_residuals}. We run 3 alternating iterations for all examples. The accumulated solve times over 3 alternating iterations for CEN and WBD optimizations are reported separately. The residuals for CoM position ($m$), EE position ($m$), and AM ($kg \cdot m^2/s$) are defined as the $l2$-norm of the error between CEN and WBD accumulated along the entire time horizon. Note that the residual for EE is also accumulated over all 4 legs. We show both the residuals from our alternating update scheme and a hierarchical one, i.e., running only one alternating iteration. It is observed that the residuals with more alternating iterations are significantly smaller and consensus performance is improved, which can be attributed to the benefit that one level of optimization iteratively receives feedback from the other.

\begin{table}[t]
    \smallskip
    \smallskip
    \caption{Solve times and residuals for trajectory generation}
    \centering
    \begin{tabular}{c c c c}
         \hline
         Desired & Solve times & Residuals for & Residuals for\\
         motions & CEN,WBD & CoM,EE,AM & CoM,EE,AM\\
          & (s) & ($m$, $kg \cdot m^2/s$) & (Hierarchical)\\
         \hline
         Cantering & 36.8, 4.9 & 1.3, 6.9, 9.9 & 1.4, 8.3, 12.7\\
         ``Parkour" & 213.0, 10.3 & 5.2, 11.3, 42.9 & 53.8, 49.7, 197.2 \\
         Inertia shaping & 35.9, 48.1 & 0, 12.9, 0 & 0, 290.0, 0\\
         Twisting jump & 23.2, 5.8 & 0.6, 1.8, 3.2 & 1.5, 2.1, 23.7\\
         \hline
    \end{tabular}
    \label{tab:solving_time_and_residuals}
\end{table}

Our offline TO pipeline is able to replicate the highly dynamic cantering gait. Motion capture (mocap) data from a real animal is used to first generate the initial seed trajectory for our robot's morphology through re-targeting \cite{RoboImitationPeng20}. Then our TO pipeline uses the re-targeted inputs to produce base and joint trajectories which are not only faithful to the original mocap recording, but also physically-correct. We benchmark our solution quality against that of only running WBD optimization. Fig.~\ref{fig:benchmark_to_wbd} shows a comparison of the front left (FL) and front right (FR) z-direction foot position trajectories for the cantering example. It is clear that WBD optimization by itself has difficulty constraining foot positions on the ground during stance phases, but our alternating pipeline can correctly enforce terrain constraints, which improves the overall solution quality.

A ``parkour" scenario, where the robot traverses over a large gap via inclined surfaces, using hopping gaits, is also successfully computed through our TO pipeline. We set up the CEN optimization by manually providing the desired contact sequence, timings, and a terrain map with a 2.4 s time horizon. Intermediate states, for each landing configuration on the inclined surfaces, are also specified as hard constraints to improve solver convergence rate. The cantering and ``parkour" examples are visualized in Fig.\ref{fig:sim_ex}.

\begin{figure}
    \centering
    \includegraphics[width=3.0in]{benchmark_to_wbd_v2-1.png}
    \caption{EE in z-axis for FL and RR legs with a flat terrain (grey), generated by using only WBD optimization (red) and our alternating one (green).}
    \label{fig:benchmark_to_wbd}
    \vspace{-0.2in}
\end{figure}

\begin{figure}[htb!]
    \centering

    \includegraphics[width=\linewidth]{zero_g_4_marked.pdf}
    \caption{Simulated motion of A1 turning $90^{\circ}$ in a zero-gravity environment }
    \label{fig:zero_g_turn}
\end{figure}

To demonstrate the benefit of momentum-inertia-aware CEN optimization in the offline TO, we study the zero-gravity body reorientation problem. Fig. \ref{fig:zero_g_turn} shows the simulation result using Drake \cite{drake}. Starting from stationary initial configuration, the robot is tasked to turn its trunk $90^{\circ}$ around the world z-axis. Since no gravity nor external contacts are present in the problem setup, the robot is forced to rely on its limb motions to produce the turning behavior. In the absence of external forces, CEN optimization by itself and all other SRBM-based approaches cannot discover meaningful EE trajectories for this specific example. However, with inertia feedback from the WBD level, the two optimization levels indeed reach dynamics consensus, as shown in Table \ref{tab:solving_time_and_residuals} with the EE residual decreasing from 290 to 12.9. This further validates that inclusion of the momentum inertia parameterization at the CEN level is physically meaningful. It is also worth noting that for this example, the WBD optimization takes significantly longer than the other examples, although the total motion horizon is short. This is because during the first alternating iteration, WBD optimization lacks good tracking reference from the CEN level, due to the absence of any contact related constraints. This further illustrates that a decent tracking reference generated from CEN level helps drastically reduce the solve time for the WBD optimization.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\linewidth]{consensus_jump_turn-1.png}
    \caption{Dynamics consensus of angular momentum for 90$^{\circ}$ twisting jump.}
    \label{fig:consensus}
\end{figure}

The $90^{\circ}$ twisting jump maneuver is also studied, shown in Fig. \ref{fig:twistingJump}. This maneuver consists of three contact phases: takeoff, aerial, and landing. The robot must gain enough AM during takeoff to achieve the goal orientation after landing. In addition to the solve times and residuals provided in Table \ref{tab:solving_time_and_residuals}, Fig. \ref{fig:consensus} illustrates desirable consensus performance, for centroidal angular momentum, between the two optimization layers. In the subsequent section, we will discuss the hardware implementation of this maneuver using the convex MPC controller proposed in Sec. \ref{sec:fullCenConvexMPC}.

\subsection{Hardware Demonstration} \label{sec:hardwareResults}
\begin{figure}[t!]
    \smallskip
    \smallskip
    \centering
    \includegraphics[width=\linewidth]{Cheetah_Jump_Turn.pdf}

    \caption{Demonstration of the Mini Cheetah executing a $90^{\circ}$ twisting jump.}
    \label{fig:twistingJump}
    \vspace{-0.2in}
\end{figure}
We validate the tracking performance of the proposed MPC controller, by reproducing the twisting jump maneuver on the Mini-Cheetah platform. Since our MPC lacks the ability to directly optimize contact timings, a heuristics-based contact detection scheme is devised to deal with the early touchdown scenario. For each leg, touchdown is estimated by detecting an abrupt velocity change for the knee joints after the takeoff phase. For each leg, if a touchdown is early, tracking reference is immediately shifted to the point of touchdown on the offline trajectory. When all four legs land, base and momentum trajectories are also shifted to their respective points of touchdown.

Our controller is deployed directly on the Mini-Cheetah onboard computer, and updates at $100$ Hz with a prediction horizon of $50$ ms. Compare to the original convex MPC controller \cite{di_carlo_dynamic_2018}, which updates at around $30$ Hz, our version updates at a higher frequency, but with a shorter prediction horizon (up to $0.5$ s prediction horizon for the original). Fig.~\ref{fig:tracking_performance_twisting_jump} demonstrates the tracking performance for the twisting jump maneuver. During the takeoff and aerial phases, both joint and momentum are tracked properly. When early touchdown is detected, tracking references are shifted heuristically as described in the previous paragraph. However, the momentum tracking during the landing phase still leaves much to be desired. As part of our future work, we expect that by properly optimizing contact timing online, instead of employing simple touchdown detection heuristics, the landing phase momentum tracking performance can be improved significantly. Nonetheless, our proposed controller successfully replicates the $90^{\circ}$ twisting jump on hardware, as shown in Fig. \ref{fig:twistingJump}.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=2.8in]{joint_AM_tracking_v4-1.png}
    \caption{Tracking performance for twisting jump by comparing the offline reference trajectories (green) and the recorded data during online execution (pink). The brown vertical dashed line denotes the early touchdown.}
    \label{fig:tracking_performance_twisting_jump}
    \vspace{-0.2in}
\end{figure}

\section{Conclusions and Future Work}
We presented a trajectory optimization and QP-based MPC framework for versatile and agile quadrupedal locomotion. The novelty of our proposed framework lied in the incorporation of both CEN and WBD during the trajectory optimization process. The offline trajectory planner adopted an alternating scheme to achieve dynamics consensus between centroidal and whole-body optimizations. The MPC controller used the generated trajectory and re-computed contact forces and joint accelerations for online execution. Our proposed pipeline has been verified through a range of dynamic motions in both simulation and real world. Future works include optimizing over foot placements and contact timings in our MPC controller to realize cantering and ``parkour" on hardware.

\let\secfnt\undefined
\newfont{\secfnt}{ptmb8t at 10pt}

\bibliographystyle{./IEEEtran}


text:
\section{Introduction}
Radio AGN are a cornerstone of modern galaxy formation models, the
current consensus favouring a paradigm in which AGN feedback
regulates star formation and shapes fundamental features of the
galaxy population like their distribution of colours and
luminosities \citep[e.g.][]{Hatton2003,Croton2006,Somerville2012}.
Observational work indicates that radio AGN, at least in the local
Universe, are more common in early-type galaxies, the detection
rate in spiral or late-type galaxies (LTGs) being around an order
of magnitude lower \citep{Ledlow2001,Kaviraj2014a}. Radio AGN are,
therefore, comparatively rare in LTGs compared to their early-type
counterparts, with relatively few examples having been studied in
the past literature
\citep[e.g.][]{Ekers1978,Ulvestad1984,Edelson1987,Ledlow1998,Morganti2001,Keel2006,Zhou2007,Yuan2008,Gliozzi2010,Inskip2010,Norris2012,Mao2014}.

\begin{figure*}
\begin{minipage}{172mm}
\begin{center}
$\begin{array}{c}
\includegraphics[width=\textwidth]{spirals.ps}
\end{array}$
\caption{SDSS colour images of a sample of VLBI-detected spiral
galaxies that underpin this study.} \label{fig:spirals}
\end{center}
\end{minipage}
\end{figure*}

The study of radio AGN in LTGs is desirable not just because they
are rare, but also because they may hold clues to the behaviour of
radio jets in disc environments. While this is uncommon at low
redshift, the bulk of the stellar and black-hole mass in today's
Universe was created at $z\sim2$
\citep[e.g.][]{Madau1998,Hopkins2006}, an epoch at which both star
formation \citep[e.g.][]{Kaviraj2013a} and black-hole growth
\cite[e.g.][]{Kocevski2012,Schawinski2012} were predominantly
hosted by late-type galaxies. The connection between the black
hole and a disc-like host system was therefore common around the
epoch of peak cosmic star formation, making radio AGN in nearby
LTGs useful laboratories for exploring this connection.

Given their relative rarity, a statistical study of radio AGN in
LTGs requires a combination of survey-scale radio data (for the
detection of large numbers of radio AGN) and high-resolution
optical imaging (for accurate determination of galaxy
morphologies). In this study, we combine mJIVE-20 (Deller \&
Middelberg 2014) with optical imaging and spectro-photometry from
the SDSS (Azabajian et al. 2009), to study radio AGN in the nearby
($z<0.3$) LTG population.

The high resolution of VLBI enables unambiguous identification of
AGN, because the $>10^6$ K temperatures required for a detection
cannot be achieved via star formation alone and requires
non-thermal sources like supernova remnants (SNRs), radio
supernovae (SNe) or AGN. Only in the very local Universe and in
vigorous starbursting systems can clusters of luminous SNe and
SNRs reach sufficiently high luminosities to be visible in a
(shallow) VLBI observation. In the extreme example of Arp 220,
which has a star formation rate of several hundred solar masses
per year \citep[e.g.][]{Iwasawa2005,Baan2007}, the brightest
VLBI-scale sources reach a peak VLBI flux density of $\sim$1
mJy/beam \citep{Lonsdale2006}. As we indicate below, our sample
consists of sources that are both more distant and have lower star
formation rates than Arp 220, making it extremely unlikely that
the relatively shallow mJIVE-20 observations (with a detection
limit of $\sim$1 mJy) would detect anything other than an AGN in
our target sample.

While in lower-resolution surveys like FIRST and NVSS, AGN
identification requires the detection of excess radio flux beyond
what is expected from star formation (biasing samples towards AGN
that dominate the star formation), or a dense gas environment
which the jet works against to produce detectable radio lobes
(introducing a bias against AGN in LTGs, which do not typically
host extensive hot gas halos), VLBI is capable of identifying AGN
irrespective of the host galaxy properties (e.g. star formation
rate, environment etc.). Some limitations of VLBI are worth noting
here. Of the total radio emission generated by AGN activity, only
a fraction will typically be confined to the parsec-scale core/jet
at the site of the central black hole, as extended radio lobes
will often be present -- hence, the compact flux fraction can be
influenced by the surrounding galactic environment. The prominence
of the compact core depends on the source age and orientation
(which can result in Doppler boosting or deboosting of the compact
emission). Finally, hotspots at the site of the jet interaction
with the interstellar medium \citep[as are seen in compact
symmetric objects -- CSOs][]{Phillips1982,Wilkinson1994} may also
be compact enough to be visible in VLBI observations, so that,
although VLBI detections can unambiguously be associated with AGN,
they cannot in every case be associated with AGN {\em cores}.
However, this last case can usually be distinguished on the basis
of morphology. Nevertheless, the unambiguity in the identification
of AGN activity makes VLBI data from mJIVE-20 the ideal route for
studying the local AGN population that inhabits LTGs.

Here, we identify our spiral AGN hosts via direct visual
inspection of the SDSS colour images of mJIVE-20 detections. We
then study their physical properties, derived from SDSS
spectro-photometric data, to explore the conditions that make it
likely for LTGs to host AGN. This is the first study where
survey-scale VLBI data is combined with visual inspection of
galaxy images to study a \emph{morphologically-selected} sample of
LTGs that host radio AGN. The sections below are organized as
follows. In Section 2, we briefly describe the mJIVE-20 and SDSS
data that underpin this study. In Section 3, we compare the
properties of LTGs that host radio AGN with those of the general
LTG population. We discuss and summarize our findings in Section
4. Throughout, we employ the WMAP7 cosmological parameters
\citep{Komatsu2011} and photometry in the AB system
\citep{Oke1983}.

\section{Data}

\subsection{Radio VLBI data}
mJIVE-20 is an ongoing survey using the VLBA that is
systematically observing objects detected by the FIRST radio
survey. mJIVE-20 uses short segments scheduled in bad weather or
with a reduced number of antennas during which no highly rated
VLBA science projects can be scheduled. The survey has targeted
$\sim$25,000 FIRST sources to date, with $\sim$5000 VLBI
detections. While the sensitivity and resolution of mJIVE-20
varies between different fields, the median detection threshold is
1.2 mJy/beam and the typical beam size is 6$\times$17
milliarcseconds. This corresponds to a detection sensitivity of
$\sim10^7$ K (the variation between fields is around a factor of
2). We refer readers to \citet{Deller2014} for further details of
the survey.

\subsection{SDSS data}
The mJIVE-20 targets are cross-matched with the latest data
release of the Sloan Digital Sky Survey
\citep[SDSS;][]{Abazajian2009}. Following \citet{Shabala2008}, we
use a matching radius of 2 arcseconds for the radio-optical
matching, which yields high (96\%) completeness and low (0.3\%)
contamination. Since morphological classification and the
identification of morphological disturbances will be an important
part of our study, we restrict ourselves to SDSS galaxies that
have photometric or spectroscopic redshifts less than 0.3. This
yields 437 SDSS galaxies with mJIVE-20 detections in this redshift
range. The entire sample is visually classified, via the SDSS
colour images, to separate early-type galaxies from LTGs. For each
galaxy, we also note the presence of morphological disturbances,
indicating that the galaxy has had a recent merger or interaction.

This visual classification yields 29 VLBI-detected LTGs. Figure
\ref{fig:spirals} presents examples of these systems. Note that
the VLBI-detected LTGs span the full spectrum of spiral
morphologies seen on the Hubble sequence
\cite[e.g.][]{Hubble1926,dev1959}, from systems that have
reasonably prominent bulges to those that are clearly dominated by
discs. It is also worth noting that the VLBI-detected LTGs show a
high incidence of morphological disturbances, suggesting that
these systems have undergone recent mergers (we return to this
point in our analysis below).

Out of the 29 VLBI-detected LTGs, 11 have SDSS spectra and
therefore spectroscopic redshifts. Our analysis below is
restricted to these objects, since quantities like absolute
magnitudes, stellar masses and local environments require an
accurate measurement of redshift. For this
spectroscopically-detected subsample, magnitudes are K-corrected
using the \texttt{KCORRECT} code of \citet{Blanton2007} and
published stellar masses, star formation rates and emission-line
classes are extracted from the latest version of the
publicly-available MPA-JHU value-added SDSS catalogue
\citep[][]{Kauffmann2003,Brinchmann2004,Tremonti2004}\footnote{http://www.mpa-garching.mpg.de/SDSS/DR7/}.

The emission-line class of each galaxy has been calculated through
a standard line-ratio analysis \citep[][see also Baldwin et al.
1981, Veilleux et al. 1987, Kewley et al. 2006]{Kauffmann2003},
using the measured values of [NII]/H$\alpha$ and [OIII]/H$\beta$.
Objects in which all four emission lines are detected with a
signal-to-noise (S/N) ratio greater than 3 are classified in the
MPA-JHU catalog as either `star-forming', `low S/N star-forming',
`composite', `Seyfert' or `LINER', depending on their location in
the [NII]/H$\alpha$ vs. [OIII]/H$\beta$ diagram. Galaxies without
a detection in all four lines are classified as `quiescent'
\citep{Kauffmann2003}.

Finally, the local environments of these
spectroscopically-detected galaxies are extracted from the group
catalogue of \citet{Yang2007}, who use a halo-based group finder
to separate the SDSS into over 300,000 structures
{\color{black}with a broad dynamic range}, {\color{black}from}
rich clusters to isolated galaxies. \citet{Yang2007} estimate the
host dark matter (DM) halo masses of {\color{black}individual}
SDSS galaxies, which are related to the traditional
classifications of environment (`field', {\color{black}`group'}
and {\color{black}`cluster'}). Cluster-sized haloes typically have
masses greater than $10^{14}$ M$_{\odot}$, while group-sized
haloes have masses between $10^{13}$ and $10^{14}$ M$_{\odot}$
\citep{Binney1987}.

\section{Radio AGN in spiral galaxies}
Given that VLBI detections are possible in rare,
intensely-starbursting galaxies like Arp 220, we first check
whether our spiral AGN hosts are consistent with an Arp 220-like
starburst. The VLBI-detected LTGs have SFRs (Figure 5) that are
significantly lower than that in Arp 220, which has an SFR of
several hundred solar masses per year
\citep[e.g.][]{Iwasawa2005,Baan2007}. They also lie much further
away than Arp 220, which has a redshift of 0.018 (Figure 3). The
VLBI radio fluxes in these systems are, therefore, inconsistent
with being driven by star formation in an Arp 220-like system.
This conclusion is consistent with the SDSS images of these
galaxies (see Figure 1), which indicate that most of these objects
are not dusty major-merger remnants like Arp 220, but are normal
massive galaxies with a disc-like morphology.

\begin{figure}
$\begin{array}{c}
\includegraphics[width=3.5in]{mjive_luminosities.ps}
\end{array}$
\caption{mJIVE-20 radio luminosities of VLBI-detected spirals
compared to that of their early-type counterparts in the same
stellar mass and redshift range.}
\label{fig:mJIVE-20_luminosities}
\end{figure}

\subsection{Physical properties of spiral AGN hosts}
We begin our analysis by studying the physical properties of the
LTGs that host radio AGN. Figure \ref{fig:mJIVE-20_luminosities}
indicates that the radio AGN in our LTG sample have similar masses
and nuclear radio luminosities as their early-type counterparts.
We are, therefore, studying radio AGN in spirals that broadly
share the same physical properties as the more commonly-found
radio AGN in early-type systems in the same stellar mass and
redshift ranges.

\begin{comment}
\begin{table*}
\begin{center}
\caption{Comparison of the properties of VLBI-detected,
VLBI-undetected and control LTGs. Columns are as follows: (1)
Galaxy number fraction on the red sequence, defined as $u-r>2.2$
(2) The fraction of galaxies that are currently merging or are
post-mergers (3-8) Fraction of galaxies in each BPT classification
class (star-forming, low S/N star forming, composites, Seyferts,
LINERs and quiescent). Errors are shown in brackets.}
\begin{tabular}{lcccccccc}\hline
                            & Red  & Merging & Star-forming & Low S/N star-forming & Composite & Seyfert & LINER & Quiescent \\\hline
    VLBI-detected LTGs      & 90\% & 53\%    & 8\%          & 8\%                  & 23\%      & 30\%    & 23\%  & 8\%\\
    VLBI-undetected LTGs    & 48\% & 28\%    & 14\%         & 7\%                  & 41\%      & 29\%    & 5\%   & 2\%\\
    Control LTGs            & 62\% & 16\%    & 12\%         & 41\%                 & 12\%      & 6\%     & 19\%  & 10\%\\

\end{tabular}
\end{center}
\label{tab:merger_fractions}
\end{table*}
\end{comment}

\begin{table*}
\begin{center}
\caption{Comparison of the properties of VLBI-detected,
VLBI-undetected and control LTGs. Columns are as follows: (1)
Galaxy number fraction on the red sequence, defined as $u-r>2.2$
(2) The fraction of galaxies that are currently merging or are
post-mergers (3-8) Fraction of galaxies in each BPT classification
class (star-forming [SF]), low S/N star forming, composites [Cp],
Seyferts [Sy], LINERs [LI] and quiescent [Qs]).}
\begin{tabular}{ccccccccc}\hline
                         & Red             & Mgr             & SF              & Low S/N SF      & Cp              & Sy              & LI              & Qs \\\hline
    VLBI-detected LTGs   & 0.90$^{(0.35)}$ & 0.53$^{(0.21)}$ & 0.08$^{(0.07)}$ & 0.08$^{(0.07)}$ & 0.23$^{(0.17)}$ & 0.30$^{(0.18)}$ & 0.23$^{(0.17)}$ & 0.08$^{(0.07)}$\\
    VLBI-undetected LTGs & 0.48$^{(0.16)}$ & 0.28$^{(0.06)}$ & 0.14$^{(0.04)}$ & 0.07$^{(0.03)}$ & 0.41$^{(0.08)}$ & 0.29$^{(0.06)}$ & 0.05$^{(0.02)}$ & 0.02$^{(0.01)}$\\
    Control LTGs         & 0.62$^{(0.06)}$ & 0.16$^{(0.02)}$ & 0.12$^{(0.02)}$ & 0.41$^{(0.04)}$ & 0.12$^{(0.02)}$ & 0.06$^{(0.02)}$ & 0.19$^{(0.03)}$ & 0.10$^{(0.02)}$\\

\end{tabular}
\end{center}
\label{tab:merger_fractions}
\end{table*}

\begin{figure}
$\begin{array}{c}
\includegraphics[width=3.5in]{mass_redshift.ps}\\
\end{array}$
\caption{Stellar mass vs redshift of VLBI-detected (red) and
VLBI-undetected (blue) and control (green) spiral galaxies. The
shaded region indicates the stellar mass range considered in our
analysis in Section 3.} \label{fig:mass_redshift}
\end{figure}

\begin{figure}
$\begin{array}{c}
\includegraphics[width=3.5in]{halomass.ps}\\
\end{array}$
\caption{Halo mass vs group-centric radii of VLBI-detected (red),
VLBI-undetected (blue) and control (green) spiral galaxies.}
\label{fig:environment}
\end{figure}

\begin{figure}
$\begin{array}{c}
\includegraphics[width=3.5in]{ur.ps}\\
\includegraphics[width=3.5in]{sfr.ps}
\end{array}$
\caption{$(u-r)$ (top) and SFR (bottom) vs stellar mass for the
VLBI-detected (red), VLBI-undetected (blue) and control (green)
spiral galaxies. The dotted line in the top panel indicates the
demarcation between the red sequence and the blue cloud
\citep[e.g.][]{Strateva2001}.} \label{fig:sfr}
\end{figure}

In our subsequent analysis, we compare the physical properties of
our VLBI-detected LTGs to two different populations of LTGs. Since
most of our VLBI-detected LTGs have stellar masses greater than
M$_* \sim 10^{10.8} M_{\odot}$ (see Figure
\ref{fig:mass_redshift}) and have been restricted to $z<0.3$, we
restrict all other samples to these stellar mass and redshift
ranges. The first population is the sample of 77 LTGs that are
mJIVE-20 targets but are undetected (we refer to these galaxies as
`VLBI-undetected LTGs'). Since the VLBI-undetected LTGs are, by
construction, detected by FIRST, we also define a `control sample'
of LTGs that are not radio detected. To construct this control
sample we select a random $\sim$230 spiral galaxies from the SDSS
via an identical visual inspection as was performed on the
mJIVE-20 targets. Figure \ref{fig:mass_redshift} presents the
stellar mass and redshift properties of each of these LTG
populations. It is worth noting that the VLBI-detected spirals
generally lie towards the upper envelope of the mass-redshift
space. KS-tests between the mass distributions of the
VLBI-detected LTGs compared to the VLBI-undetected and control
LTGs yield p-values of $\sim$10$^{-5}$ and $\sim$10$^{-4}$
respectively, indicating that the spirals with AGN are likely to
be drawn from different parent distributions. The offset towards
higher masses in the VLBI-detected systems is plausibly due to
larger galaxies having bigger black holes \citep[][]{Gultekin2009}
which, in turn, produce more powerful jets
\citep[e.g.][]{Shabala2008}.

\begin{figure}
$\begin{array}{c}
\includegraphics[width=3.5in]{bpt_vlbidetectedspirals.ps}\\
\includegraphics[width=3.5in]{bpt_vlbiundetectedspirals.ps}\\
\includegraphics[width=3.5in]{bpt_controlspirals.ps}
\end{array}$
\caption{BPT classifications of VLBI-detected (top),
VLBI-undetected (middle) and control (bottom) spiral galaxies. The
fractions of galaxies in each morphological class are indicated in
Table 1.} \label{fig:bpt}
\end{figure}

Figure \ref{fig:environment} shows the environments of the spiral
AGN hosts compared to that of their VLBI-undetected and control
counterparts. Paired (i.e. 2D) KS-tests in the halo mass - group
centric radius plane between the VLBI-detected LTGs and their
VLBI-undetected and control counterparts yield p-values of 0.06
and 0.01 respectively. Thus, while the environments of the
VLBI-detected and undetected populations are likely to be similar,
there is a hint that the VLBI-detected LTGs are somewhat different
from the control LTGs. The difference is probably driven by the
fact that the VLBI detections appear to avoid the the highest halo
masses compared to the control sample. However, local environment,
on the whole, appears not to play a major role in determining the
presence of radio AGN in LTGs.

The VLBI-detected spirals generally show redder $(u-r)$ colours
and correspondingly lower SFRs than the undetected and control
populations (Figure \ref{fig:sfr}). KS-tests between the
VLBI-detected and control LTGs in either of these quantities yield
very low p-values ($<$10$^{-3}$), indicating that the spiral AGN
hosts are likely to be different from the general LTG population
in these properties. $\sim$90\% of the VLBI-detected LTGs are on
the red sequence, compared to $\sim$50\% of the VLBI-detected
sample and $\sim$60\% of the control LTG sample (column 1 in Table
1). In a similar vein, only 2 out of the 9 VLBI-detected LTGs lie
on the SFR-mass parameter space defined by the bulk of the control
and VLBI-undetected populations \footnote{It is worth noting here
that some of our spiral AGN are on the star formation main
sequence, illustrating the utility of VLBI in being able to detect
AGN regardless of the level of star formation activity in the host
system.}.

As has been noted in the literature \citep[e.g.][]{Greene2013},
the mechanism that dissipates angular momentum and allows gas to
accrete on to the central black hole remains unclear. The removal
of angular momentum could take place as a result of processes such
as merging \citep[e.g.][]{Sanders1988}, or via circumnuclear
structures such as bars \citep[e.g.][]{Maciejewski2002,Kim2012} or
nuclear spirals
\citep[e.g.][]{Englmaier2000,Martini2003,Davies2009}. Column 2 in
Table 2 presents merger fractions in our three LTG samples. The
VLBI-detected LTGs show the highest merger fraction ($\sim$50\%),
which is around a factor of 2 higher than that in their undetected
counterparts and around a factor of 4 higher than that in the
control LTGs. While other processes that drive the dissipation of
angular momentum cannot be completely ruled out, the remarkably
high merger fraction in the VLBI-detected spirals suggests that
mergers are likely to play a significant role in triggering radio
AGN in spiral hosts. It is worth noting that the SDSS images that
form the basis of the visual inspection and the identification of
morphological disturbances are relatively shallow, with standard
exposure times of 54 seconds. As noted in the recent literature
\citep[e.g.][]{Kaviraj2010}, merger fractions derived from such
images are, therefore, strictly lower limits. It is possible (and
likely) that a much larger fraction of the VLBI-detected LTGs
carry morphological disturbances, that are invisible in the
shallow standard-depth SDSS images.

A merger trigger for the AGN appears to be supported by the
emission-line analysis of our sample (Figure \ref{fig:bpt} and
columns 4-8 in Table 2). The BPT analysis suggests an increasing
fraction of optical emission-line AGN as we transition from the
control LTG sample to the VLBI-detected population. While
$\sim$37\% of the control sample show evidence for an optical AGN
(i.e. systems classified as composite, Seyferts or LINERs), the
corresponding fraction in the VLBI-detected population is
$\sim$76\%. It is worth noting that in early-type galaxies, some
LINERs may not be driven by AGN \citep[e.g.][]{Sarzi2010}. Even if
this was also the case in LTGs and the LINERs were removed from
this argument, the VLBI-detected spirals still show a surfeit of
optical AGN of around a factor of 3. This high fraction of optical
AGN appears consistent with a radiatively-efficient cold-mode type
accretion, which might be expected if the gas has been accreted
via a merger, as seems to be the case for our sample
\citep[e.g.][]{Hardcastle2007,Best2012,Shabala2012}. Our results
appear consistent with recent work that has suggested a strong
connection between merging and the onset of AGN activity at low
redshift
\citep[e.g.][]{Koss2010,Ellison2011,Scott2014,Kaviraj2014a},
although the role of mergers in triggering AGN may become
negligible at high redshift
\citep[e.g.][]{Chiaberge2011,Schawinski2012,Kocevski2012}.

\begin{figure*}
\begin{minipage}{172mm}
\begin{center}
$\begin{array}{ccc}
\includegraphics[width=2.1in]{5850.dr10.jpg.ps}  & \includegraphics[width=2.1in]{5850_first.bmp.ps}
\includegraphics[width=2.1in]{MJV16427.clean2.ps}\\
\includegraphics[width=2.1in]{22008.dr10.jpg.ps}  & \includegraphics[width=2.1in]{22008_first.bmp.ps}
\includegraphics[width=2.1in]{MJV24477.clean2.ps}\\
\includegraphics[width=2.1in]{8194.dr10.jpg.ps}  & \includegraphics[width=2.1in]{8194_first.bmp.ps}
\includegraphics[width=2.1in]{MJV01925.clean2.ps}
\end{array}$
\caption{SDSS (left), FIRST (middle) and mJIVE-20 (right) images
of three spiral AGN hosts in our sample, that are close to edge-on
with evidence for a jet in the VLBI image. Note that, in the last
example (bottom row), the jet is shown slightly off centre, in
order to capture all of the structure around this galaxy on the
same spatial scales as for the other two
galaxies.}
\label{fig:VLBI_jets}
\end{center}
\end{minipage}
\end{figure*}

The analysis above indicates that local LTGs that host radio AGN
(with luminosities similar to those found in their early-type
counterparts) tend to be massive spirals with low SFRs, with
mergers being a plausible trigger for the onset of the AGN. It is
worth noting that that the requirements for LTGs to host radio AGN
are akin to those in their early-type counterparts that inhabit
similar (low-density) environments i.e. massive galaxies that have
red colours and are likely to be involved in a merger
\citep{Kaviraj2014a}. Nevertheless, the detection rate in
early-types is almost an order of magnitude higher
\citep{Kaviraj2014a}. While a thorough investigation is beyond the
scope of this particular paper, we speculate on some of the
potential reasons for this discrepancy. It is possible that hot
gas reservoirs in early-type galaxies could provide an extra
source of fuel for the black hole - we note, however, that the
\citet{Kaviraj2014a} sample does not contain central cluster
galaxies where hot gas cooling is expected to play a significant
role in AGN triggering
\citep[e.g.][]{Tabor1993,McNamara2007,Cattaneo2009,Fabian2012}.
Thus, it is unclear if the observed discrepancy in detection rate
can be completely explained solely via this argument.
Alternatively (or in addition), black hole masses may trace bulge
mass more closely than total galaxy stellar mass, implying that
even massive LTGs do not typically have black holes as massive as
those in their early-type counterparts. Finally, the `spin
paradigm', which postulates that radio jets are launched by
rapidly spinning black holes
\citep[e.g.][]{Blandford1977,Wilson1995,Sikora2007,Bagchi2014}
that are expected to reside primarily in bulge-dominated galaxies,
may also play a leading role in explaining the surfeit of radio
AGN in early-type systems compared to their late-type
counterparts.

\subsection{Implications for AGN feedback in spiral AGN hosts}
The recent literature indicates that, in low-density environments
in the \emph{nearby} Universe, AGN overwhelmingly avoid systems
that are in the blue cloud i.e. that are experiencing the peak of
their star-formation episode
\citep[e.g.][]{Schawinski2007,Wild2010,Carpineti2012,Kaviraj2011,Kaviraj2014b}.
Radio AGN in early-types, for example, cluster strongly on the red
sequence, where star formation rates have declined significantly
\citep{Kaviraj2014a}. Given that typical AGN lifetimes \citep[a
few times 10$^7$ yr, e.g.][]{Shabala2008} are much shorter than
the transit time from blue cloud to red sequence, this suggests a
time lag between the onset of star formation and the triggering of
AGN. Interestingly, this time lag shows no dependence on the
intensity of star formation (and therefore the gas fraction) in
the host system and is observed across the entire spectrum of star
formation activity, from intensely-starbursting luminous infrared
galaxies \citep{Kaviraj2009a} to the relatively weakly
star-forming early-type systems
\citep[e.g.][]{Wild2010,Shabala2012}. The implication of this time
delay is that the AGN typically does not couple to the host's gas
reservoir during the main starburst phase, thus reducing its
ability to strongly regulate the stellar mass growth.

Our study reinforces this picture, by extending these results to
spiral AGN hosts. In a similar vein to what is found in early-type
galaxies, radio AGN in the nearby LTG population preferentially
inhabit systems on the red sequence that have low SFRs. Unless the
typical AGN lifetime (i.e. the duration of the `on' phase) is
significantly longer in LTGs compared to their early-type
counterparts, the radio AGN in our spiral hosts are also likely to
have been triggered after the peak of the star formation activity
has elapsed in these systems, limiting their ability to impart
feedback on the gas reservoir. Taken together with the recent
literature, our results indicate that the time delay between the
onset of merger-driven star formation and the subsequent
triggering of AGN is likely to be a ubiquitous feature of radio
AGN in the nearby Universe.

\subsection{Jet orientations in spiral AGN hosts}
We conclude this section by exploring the morphology of the jets
in our VLBI-detected spirals. Three galaxies in our sample show
clear VLBI-scale jets. In Figure \ref{fig:VLBI_jets}, we present
the VLBI radio maps of these galaxies, together with their optical
SDSS and arcsecond-scale radio images from FIRST. In all cases,
the VLBI-scale jets are clearly close to being orthogonal to the
galaxy disc. Since the jet is expected to align with the spin axis
of the black hole \citep[e.g.][]{McKinney2013}, the orthogonality
observed in our galaxies is not unexpected, since co-evolution of
the black hole and the host spiral should produce alignment
between the black hole's spin axis and the stellar disc. In such a
scenario one would expect the jet to be orthogonal to the major
axis of the spiral galaxy, as is observed in our galaxies.

It is worth recalling that a high fraction (50\%+) of our LTG AGN
are likely merger remnants. One might expect that galaxy mergers
could disrupt the alignment between the black hole spin axis and
the stellar disc e.g. due to the coalescence of two similar mass
black holes in a major merger. However, it is important to note
that the interactions in question here are, by definition,
\emph{minor} mergers (i.e. mergers with small satellites), because
the discs in our LTG AGN remain intact, contrary to what might be
expected in major-merger remnants
\citep[e.g.][]{Barnes1992a,Springel2005a}. Our results therefore
suggest that minor mergers are unlikely to disrupt the alignment
between the black hole spin axis and the stellar disc (which may
naturally come about due to secular co-evolution of these two
systems)\footnote{Note that the well-known misalignment of
kpc-scale radio jets with the galaxy disk \citep{Ulvestad1984}
does not pose a serious problem to this interpretation, since
parsec and kpc-scale radio jets are often spectacularly
misaligned, as illustrated in Cen A \citep{Israel1998,Feain2009}.}

Finally, it is worth noting that several authors have studied the
relative geometry of galaxy-scale and circumnuclear structures in
systems that host AGN. Most recently, \citet{Greene2013} used
H$_2$O megamaser observations to map the circumnuclear discs in
seven nearby Seyfert galaxies, and found that while the maser
discs are oriented perpendicular to the direction of VLBI-scale
radio jets (consistent with jet generation models in which the
magnetic field associated with the black hole's accretion disk
determines jet orientation, e.g. \citet{McKinney2013}), the maser
and galaxy discs did not align in most of their sample. While the
reasons for this discrepancy with our results are difficult to
probe further using the available data, several causes for such
misalignments are possible, including disk warping on parsec
scales, subparsec-scale torques from stars in a cusp around the
black hole and natural changes in the angular momentum vector of
(secularly) infalling gas as it is torqued by the stars
\citep[e.g.][]{Maloney1996,Pringle1997,Hopkins2010,Gammie2000,Greene2013}.

\section{Summary and discussion}
We have combined the mJIVE-20 radio VLBI survey with the SDSS to
study the rare population of radio AGN in LTGs that have radio
luminosities similar to the more commonly-found early-type AGN in
the same stellar mass and redshift ranges. As noted in the
introduction, the study of radio AGN in LTGs is desirable, not
just because they are rare, but also because they are good
laboratories for exploring the behaviour of radio jets in disc
environments. While this is relatively rare at low redshift, the
connection between black holes and discs is likely to have been a
common feature of galaxy evolution around the peak epoch of star
formation ($z\sim2$) when most of today's stellar mass was
created. Radio AGN in LTGs may therefore hold clues to the
interplay between AGN and star formation at these critical epochs.
Our main results are as follows:

\begin{itemize}

    \item The radio AGN in LTGs are predominantly found in galaxies that have high stellar masses (M$_* > 10^{10.8}
    M_{\odot}$), red colours and low star formation rates, with a negligible dependence on the detailed morphology or local environment of the host
    galaxy.

    \item  A high fraction ($\sim50$\%) of spiral AGN hosts
    exhibit morphological disturbances that are indicative of recent mergers. The
    merger fraction is around a factor of 4 higher than that in
    the general LTG population, indicating that mergers are likely to be an
    important trigger for the radio AGN in these systems.

    \item The fraction of LTG AGN hosts that reside on the red sequence is
    almost a factor of 2 higher than that in the general LTG population. In a similar vein
    to what has been suggested in the recent literature, the predominantly red colours of the LTG AGN
    hosts suggest that the triggering of radio AGN in these systems
    is generally delayed with respect to the peak of the associated star formation
    episode, reducing the AGN's ability to regulate the stellar mass growth. Together with the recent
    literature, it seems reasonable to suggest that in starbursts
    where the the gas is accreted via a merger in the nearby Universe, AGN feedback is unlikely
    to strongly regulate the resultant star formation.

    \item In LTGs that exhibit parsec-scale jets in the VLBI images,
    we find that the orientations of these jets appear to be roughly perpendicular to the
    major axis of the host galaxy. If the black holes and discs
    co-evolve in these systems, then one might expect the black
    hole spin axes to align with the stellar discs which, in turn, would produce the observed alignment between the galaxy disc
    and the jet.

\end{itemize}

\section*{Acknowledgements}
We are grateful to Martin Hardcastle for many illuminating
discussions. We thank Meg Urry and Kevin Schawinski for comments.
SK is grateful for support from the University of Tasmania (UTas)
via a UTas Visiting Scholarship and acknowledges a Senior Research
Fellowship from Worcester College Oxford. SSS acknowledges an ARC
Early-Career Fellowship (DE130101399). ATD was supported by an NWO
Veni Fellowship.

\bibliographystyle{mn2e}

redpajama_set_name: RedPajamaArXiv