11institutetext: Max Planck Institute for Intelligent Systems, Tuebingen,
Germany,
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>
22institutetext: LAAS-CNRS, Université de Toulouse, CNRS, Toulouse,
33institutetext: Industrial Engineering Department, University of Trento,
Italy,
33email<EMAIL_ADDRESS>44institutetext: Tandon School of
Engineering, New York University, New York, USA,
44email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Nonlinear Stochastic Trajectory Optimization for Centroidal Momentum Motion
Generation of Legged Robots

Ahmad Gazar 11 Majid Khadiv 11 Sébastien Kleff 2244 Andrea Del Prete 33
Ludovic Righetti 1144

###### Abstract

Generation of robust trajectories for legged robots remains a challenging task
due to the underlying nonlinear, hybrid and intrinsically unstable dynamics
which needs to be stabilized through limited contact forces. Furthermore,
disturbances arising from unmodelled contact interactions with the environment
and model mismatches can hinder the quality of the planned trajectories
leading to unsafe motions. In this work, we propose to use stochastic
trajectory optimization for generating robust centroidal momentum trajectories
to account for additive uncertainties on the model dynamics and parametric
uncertainties on contact locations. Through an alternation between the robust
centroidal and whole-body trajectory optimizations, we generate robust
momentum trajectories while being consistent with the whole-body dynamics. We
perform an extensive set of simulations subject to different uncertainties on
a quadruped robot showing that our stochastic trajectory optimization problem
reduces the amount of foot slippage for different gaits while achieving better
performance over deterministic planning.

###### keywords:

stochastic optimal control, chance-constraints, trajectory optimization,
legged robots

##  1 Introduction

Trajectory optimization has become a dominant paradigm for planning and
control of whole-body trajectories for legged robots [1, 2, 3, 4]. Despite
recent advancements in real-time trajectory re-planning [5, 6], dealing with
uncertainties and assessing risk remains an open challenge in controlling
legged robots. Uncertainties arising from inaccurate modelling, estimation
errors, etc., might cause the robot to make/break contact at a different
contact position and/or timing than the planned trajectory. This often leads
to constraint violations in the form of slipping or collision of the robot
with its environment causing safety hazards.

Figure 1: Robust trajectory optimization framework alternating between
centroidal states of whole-body DDP motions and stochastic centroidal SCP
motions.

Most of the work that took into account uncertainties in the control of legged
locomotion mainly resorts to linear models, such as the Linear Inverted
Pendulum Model (LIPM). For instance, Villa et al. [7] used a tube-based linear
Robust Model Predictive Control (RMPC) to plan Center of Mass (CoM)
trajectories subject to worst-case disturbances on the Center of Pressure
(CoP). To reduce the conservativeness of RMPC, Gazar et al. [8] resorted to
linear stochastic MPC subject to stochastic additive disturbances on the
dynamics and CoM linear chance-constraints. Yeganegi et al. [9] used Bayesian
Optimization to learn cost function weights to achieve robust walking motions
under different uncertainties. However, assuming fixed height for the CoM and
zero angular momentum around the CoM limits the range of motions and cannot
plan trajectories in non co-planar multi-contact scenarios.

There has been some effort in the literature to take into account
uncertainties in the nonlinear trajectory optimization problem for legged
locomotion. In these approaches, Differential Dynamic Programming (DDP) has
been the method of choice as it takes into account the structure of the
underlying optimal control problem to find the optimal trajectories and
control policy [10, 11]. In an effort to take uncertainties into account,
Mordatch et al. [12] used an offline ensemble DDP enabling them to transfer
the whole-body trajectories to a real humanoid robot. Recently, Hammoud et al.
[13] used a risk-sensitive version of DDP to optimize for impedance gains to
deal with uncertainties in contact. Despite the appealing convergence
properties of DDP-based approaches, inequality constraints can only be
considered implicitly in the form of penalties, which makes it difficult to
reason about the effect of uncertainties on the constraint satisfaction
robustness. Tangent to that work, Darnach et al. [14] used Expected Residual
Minimization (ERM) to solve a contact-implicit stochastic complementarity
problem. However, contact-implicit trajectory optimization can be quite hard
to tune, computationally demanding, and gets easily stuck in local minimas. A
middle-ground model between whole body and template models for legged
locomotion is the centroidal momentum dynamics [15]. While this model captures
the dynamics between the external forces and the centroidal states exactly
[16], it is relatively low dimensional and the nonlinearity structure in the
model is well understood [17]. Hence, it is more tractable to use direct
methods with explicit consideration of the force and CoM constraints and
design force trajectories that are robust to the different types of
uncertainties in the model. In this work, we propose to solve stochastic
trajectory optimization problems for centroidal momentum generation of legged
robots in the presence of additive uncertainties in the model dynamics and
uncertainties in the contact location. We write down the friction cone
constraints as chance constraints and study the effect of the uncertainties on
the satisfaction of these constraints. Our work is related to [18, 19], where
the motion of the CoM was restricted within certain bounds to make the
underlying optimal control problem convex. In our work instead, we take into
account different types of realistic uncertainties in the trajectory
optimization problem to find trajectories that are robust against them and are
better transferable to the real world. Moreover, compared to [18, 19], we
construct the problem as a stochastic optimization rather than Min-max
formulation which is a less conservative formulation and facilitates the
generation of more complex behaviours.

Our main contributions are: 1) We propose to use stochastic trajectory
optimization using Sequential Convex Programming (SCP) for generating robust
centroidal momentum trajectories subject to additive uncertainties in the
dynamics, as well as parametric contact position uncertainties. 2) By
considering chance constraints on the friction cones, this is the first work
(to the best of our knowledge) that designs controllers for legged robots that
generates robust force trajectories subject to contact location uncertainty.
In the same spirit as [17], we propose a whole-body trajectory planning
framework that alternates between whole-body motion planning and centroidal
momentum planning. Contrary to [17]—but similarly to [20]—we use whole-body
DDP instead of a kinematics optimization, and use the solution to warm-start
the stochastic SCP centroidal momentum solver. The resulting (robust) momentum
trajectories are later used by the whole-body DDP again to track the resulting
robust centroidal trajectories and contact force trajectories. This way, we
make sure that the generated momentum trajectories are robust while being
consistent with the robot whole-body dynamics. 3) Finally, we run 400 Monte-
carlo simulations on the open-source quadruped robot Solo [21] in a Pybullet
simulation environment [22] for dynamic trotting and bounding gaits while
applying different disturbances. We show that Stochastic trajectory
optimization is able to complete all the motions safely while reducing feet
slippage, and achieving better centroidal tracking performance over the
deterministic planning 111submission video..

##  2 Background

###  2.1 Notation

A random variable $x$ following a distribution $\mathcal{Q}$ is denoted as
$x\sim\mathcal{Q}$, with $\mathbb{E}[x]$ being the expected value of $x$, and
$\boldsymbol{\Sigma}_{x}\triangleq\mathbb{E}[(\boldsymbol{x}-\mathbb{E}[\boldsymbol{x}])(\boldsymbol{x}-\mathbb{E}[\boldsymbol{x}])^{\top}]$.
The weighted $l_{2}$ norm is denoted as
$\mathinner{\\!\left\lVert\boldsymbol{y}\right\rVert}_{\boldsymbol{P}}\triangleq\boldsymbol{y}^{\top}\boldsymbol{P}\boldsymbol{y}$.

###  2.2 Robot Dynamics

The full-body dynamics of a floating-base robot in contact with the
environment can be derived using Euler-Lagrange equations of motion as follows
[16]:

$\displaystyle\boldsymbol{M}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{h}(\boldsymbol{q},\dot{\boldsymbol{q}})=\sum^{n_{c}}_{i=0}\boldsymbol{J}^{\top}_{e,i}(\boldsymbol{q})\boldsymbol{\lambda}_{e,i}+\boldsymbol{S}^{\top}\boldsymbol{\tau}_{q},$
(1)

where
$\boldsymbol{q}=\begin{bmatrix}\hat{\boldsymbol{x}}^{\top},\hat{\boldsymbol{q}}^{\top},\boldsymbol{\theta}^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{3}\times\mathbb{SO}(3)\times\mathbb{R}^{n}$
represents the generalized robot position characterizing the robot’s floating
base pose (position and orientation) w.r.t. an inertial frame $\mathcal{I}$,
and the joint positions respectively.
$\boldsymbol{M}(\boldsymbol{q})\in\mathbb{R}^{(n+6)\times(n+6)}$ denotes the
inertia matrix, and
$\boldsymbol{h}(\boldsymbol{q},\dot{\boldsymbol{q}})\in\mathbb{R}^{n+6}$ is
the vector capturing the Coriolis, centrifugal, gravity and joint friction
forces. $\boldsymbol{J}_{e,i}$ is the associated jacobian of the $i$-th end-
effector wrench $\boldsymbol{\lambda}_{e,i}$ acting on the environment.
Finally, $\boldsymbol{S}=\begin{bmatrix}\boldsymbol{0}_{(n\times
6)},\,\boldsymbol{I}_{n}\end{bmatrix}$ is the selector matrix of the actuated
joint torques $\boldsymbol{\tau}_{q}$. We can further split (1) into its
under-actuated and actuated parts;

$\displaystyle\boldsymbol{M}^{u}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{h}^{u}(\boldsymbol{q},\dot{\boldsymbol{q}})$
$\displaystyle=\sum^{n_{c}}_{i=0}\boldsymbol{J}^{u^{\top}}_{e,i}(\boldsymbol{q})\boldsymbol{\lambda}_{e,i},$
(2a)
$\displaystyle\boldsymbol{M}^{a}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\boldsymbol{h}^{a}(\boldsymbol{q},\dot{\boldsymbol{q}})$
$\displaystyle=\sum^{n_{c}}_{i=0}\boldsymbol{J}^{a^{\top}}_{e,i}(\boldsymbol{q})\boldsymbol{\lambda}_{e,i}+\boldsymbol{\tau}_{q}.$
(2b)

By writing down the floating-base dynamics for the CoM instead of the floating
base position $\hat{x}$, we obtain the following relationship between the
centroidal momentum dynamics $\dot{\boldsymbol{h}}_{\mathcal{G}}$ and the
generalized velocities $\dot{\boldsymbol{q}}$

$\dot{\boldsymbol{h}}_{\mathcal{G}}=\begin{bmatrix}\dot{\boldsymbol{\kappa}}\\\
\dot{\boldsymbol{l}}\end{bmatrix}=\dot{\boldsymbol{A}}_{\mathcal{G}}(\boldsymbol{q})\ddot{\boldsymbol{q}}+\dot{\boldsymbol{A}}_{\mathcal{G}}(\boldsymbol{q})\dot{\boldsymbol{q}},$
(3)

via the Centroidal Momentum Matrix (CMM)
$\boldsymbol{A}_{\mathcal{G}}\in\mathbb{R}^{6\times(n+6)}$ [15]. The angular
and linear momenta are denoted as $\boldsymbol{\kappa}$ and
$\boldsymbol{l}\in\mathbb{R}^{3}$ respectively. Given (3), we are interested
in planning desired centroidal momentum trajectories that satisfy the
following Newton-Euler dynamics:

$\dot{\boldsymbol{h}}=\begin{bmatrix}\sum^{n_{c}}_{i=0}(\boldsymbol{p}_{e,i}+\boldsymbol{R}_{e,i}^{x,y}\boldsymbol{\zeta}_{e,i}-\boldsymbol{c})\times\boldsymbol{f}_{e,i}+\boldsymbol{R}^{z}_{e,i}\tau_{e,i}\\\
m\boldsymbol{g}+\sum^{n_{c}}_{i=0}\boldsymbol{f}_{e,i}\end{bmatrix}$ (4)

where $\boldsymbol{c}\in\mathbb{R}^{3}$ represents the robot’s CoM, at which
the total mass $m$ of the robot is concentrated.
$\boldsymbol{p}_{e,i}\in\mathbb{R}^{3}$ is the $i$-th end-effector’s contact
position, with $\boldsymbol{\zeta}_{e,i}\in\mathbb{R}^{2}$ being the local
Center of Pressure (CoP). $\boldsymbol{f}_{e,i}\in\mathbb{R}^{3}$ and
$\tau_{e,i}\in\mathbb{R}$ represent end-effector’s contact forces and torque
for flat-footed robots, respectively. The rotation matrix
$\boldsymbol{R}_{e,i}\in\mathbb{SO}(3)$ captures the contact normals mapping
quantities from the $i$-th end-effector’s frame to the inertial frame. Note
that for point-footed robots, which we consider from now on,
$\boldsymbol{\zeta}_{e,i}$ and $\tau_{e,i}$ are always null, but the same
analysis still holds for flat-footed robots.

###  2.3 Centroidal Momentum Trajectory Optimization

First we present the deterministic nonlinear discrete-time optimal control
problem (OCP) for centroidal momentum trajectory optimization with fixed
contact position and timing.

######  Problem 1.

Nominal Optimal Control Problem (NOCP)

$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{u}\end{subarray}}.\,\,\,l_{f}(\boldsymbol{x}_{N})+\sum^{N-1}_{i=0}l(\boldsymbol{x}_{k},\boldsymbol{u}_{k})$
(5a) s.t. $\displaystyle\quad\begin{bmatrix}\boldsymbol{c}_{k+1}\\\ \\\
\boldsymbol{l}_{k+1}\\\ \\\
\boldsymbol{\kappa}_{k+1}\end{bmatrix}=\begin{bmatrix}\boldsymbol{c}_{k}+\frac{1}{m}\boldsymbol{l}_{k}\Delta_{k}\\\
\\\
\boldsymbol{l}_{k}+m\boldsymbol{g}\Delta_{k}+\sum^{n_{c}}_{i=0}\boldsymbol{f}_{{e,i}_{k}}\Delta_{k}\\\
\\\
\boldsymbol{\kappa}_{k}+\sum^{n_{c}}_{i=0}(\boldsymbol{p}_{{e,i}_{k}}-\boldsymbol{c}_{k})\times\boldsymbol{f}_{{e,i}_{k}}\Delta_{k}\end{bmatrix},$
(5b)
$\displaystyle\quad-\mu\mathfrak{f}^{z}_{{e,i}_{k}}\leq\mathfrak{f}^{x}_{{e,i}_{k}}\leq\mu\mathfrak{f}^{z}_{{e,i}_{k}},\quad\mathfrak{f}^{z}_{{e,i}_{k}}\geq
0,$ (5c)
$\displaystyle\quad-\mu\mathfrak{f}^{z}_{{e,i}_{k}}\leq\mathfrak{f}^{y}_{{e,i}_{k}}\leq\mu\mathfrak{f}^{z}_{{e,i}_{k}},\quad\mathfrak{f}^{z}_{{e,i}_{k}}\geq
0,$ (5d)
$\displaystyle\quad|\boldsymbol{p}_{{e,i}_{k}}-\boldsymbol{c}_{k}|\leq\mathcal{L}^{\max}_{e,i},$
(5e) $\displaystyle\quad\boldsymbol{x}_{0}=\boldsymbol{x}(0),$ (5f)
$\displaystyle\quad\boldsymbol{x}_{f}=\boldsymbol{x}(N),$ (5g)
$\displaystyle\quad\forall k\in\\{0,1,..,N-1\\},$ (5h)

where $\boldsymbol{x}=(\boldsymbol{x}_{0},\dots,\boldsymbol{x}_{N})$ with
$\boldsymbol{x}_{k}\in\mathbb{R}^{9}=(\boldsymbol{c}_{k},\boldsymbol{l}_{k},\boldsymbol{\kappa}_{k})$,
and $\boldsymbol{u}=(\boldsymbol{u}_{0},\dots,\boldsymbol{u}_{N-1})$ with
$\boldsymbol{u}_{k}\in\mathbb{R}^{3n_{c}}=(\boldsymbol{f}_{e,0},\dots,\boldsymbol{f}_{e,n_{c}})$
are the states and control optimizers along the control horizon N. The
centroidal momentum dynamics are discretized with a time-step $\Delta_{k}$
using an explicit Euler integration scheme (5b), where (5f)-(5g) represents
the initial and final conditions respectively. In order to avoid contact
slippage, the local contact forces in the end-effector frame
($\boldsymbol{\mathfrak{f}}=\boldsymbol{R}^{T}\boldsymbol{f}$) are constrained
inside the linearized friction cone constraints (5c)-(5d), where the static
coefficient of friction is denoted as $\mu$ with the vertical component of the
force being positive. Finally, the CoM is constrained to be within the leg
length reachability limits (5e).

##  3 Stochastic Optimal Control for Centroidal Momentum Trajectory
Optimization

In this section we present a stochastic version of problem (1) that takes into
account additive stochastic uncertainties on the centroidal momentum dynamics
as well as contact position uncertainties subject to friction pyramid chance
constraints. We consider the following discrete-time stochastic nonlinear OCP:

######  Problem 2.

Stochastic Optimal Control Problem (SOCP)

$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{x},\boldsymbol{u}\end{subarray}}.\,\,\,l_{f}(\boldsymbol{x}_{N})+\sum^{N-1}_{i=0}l(\boldsymbol{x}_{k},\boldsymbol{u}_{k})$
(6a) s.t.
$\displaystyle\quad\boldsymbol{x}_{k+1}=\boldsymbol{f}(\boldsymbol{x}_{k},\boldsymbol{u}_{k},\boldsymbol{\theta}_{k},\boldsymbol{w}_{k}),$
(6b)
$\displaystyle\quad\text{Pr}(\boldsymbol{H}\boldsymbol{x_{k}}\leq\boldsymbol{h})\geq\alpha_{x},$
(6c)
$\displaystyle\quad\text{Pr}(\boldsymbol{G}\boldsymbol{u}_{k}\leq\boldsymbol{g})\geq\alpha_{u},$
(6d) $\displaystyle\quad\boldsymbol{x}_{0}=\boldsymbol{x}(0),$ (6e)
$\displaystyle\quad\boldsymbol{x}_{f}=\boldsymbol{x}(N),$ (6f)
$\displaystyle\quad\forall k\in\\{0,1,..,N-1\\}.$ (6g)

With an abuse of notation from Problem (1), $\boldsymbol{x}_{k}$ and
$\boldsymbol{u}_{k}$ will be considered the stochastic state and control
policies evolving according to the parametric and additive stochastic
disturbance realizations $\boldsymbol{\theta}_{k}$, and $\boldsymbol{w}_{k}$.
(6c)-(6d) are the state and control polytopic joint chance constraints with
$\alpha_{x}$ and $\alpha_{u}$ being the probability levels of state and
control constraint satisfaction respectively.

######  Assumption 1.

(i.i.d. Gaussian disturbances)
$\boldsymbol{\theta}_{k}\sim\mathcal{N}(\mathbb{E}[\boldsymbol{\theta}_{k}],\boldsymbol{\Sigma}_{\boldsymbol{\theta}_{k}})$,
and
$\boldsymbol{w}_{k}\sim\mathcal{N}(\mathbb{E}[\boldsymbol{w}_{k}],\boldsymbol{\Sigma}_{\boldsymbol{w}_{k}})$
are assumed to be independent and identically distributed (i.i.d.) disturbance
realisations following Gaussian distributions.
$\mathbb{E}[\boldsymbol{\theta}_{k}]=\boldsymbol{p}_{{e,i}_{k}}$, and
$\boldsymbol{\Sigma}_{\boldsymbol{\theta}_{k}}\in\mathbb{R}^{(3n_{c}\times
3n_{c})}$ represent the mean and covariance of the contact positions
respectively. $\mathbb{E}[\boldsymbol{w}_{k}]=\boldsymbol{0}$ and
$\boldsymbol{\Sigma}_{\boldsymbol{w}_{k}}\in\mathbb{R}^{9\times 9}$ are the
mean and covariance of the additive noise on the centroidal dynamics.

###  3.1 Individual Chance Constraints Reformulation

Solving the above joint chance constraints (6c)-(6d) involves the integration
of multi-dimensional Gaussian Probability Density Functions (PDFs), which
becomes computationally intractable for high dimensions. One effective
solution is to use Boole’s inequality:

$\text{Pr}(\bigvee^{n}_{i=1}\boldsymbol{C}_{i})\leq\sum^{n}_{i}\text{Pr}(\boldsymbol{C}_{i})$
(7)

as a conservative union bound on the joint chance constraints [23]. We can
rewrite the complement of the state chance constraints as a conjunction of
individual chance constraints as follows:

$\displaystyle(\ref{eq:state chance constraints})$
$\displaystyle=\text{Pr}(\bigwedge^{l_{x}}_{i=1}\boldsymbol{H}_{i}\boldsymbol{x}\leq
h_{i})\geq\alpha_{x},$

which can be written conservatively as

$\displaystyle\text{Pr}(\bigvee^{l_{x}}_{i=1}\boldsymbol{H}_{i}\boldsymbol{x}>h_{i})\leq
1-\alpha_{x}.$ (8)

By applying Boole’s inequality on the above equation, and allocating
constraint violation risk equally $\epsilon_{x_{i}}=(1-\alpha_{x})/l_{x}$,
with $l_{x}$ being the number of intersecting hyper-planes forming the state
joint polytopic constraint, we reach

$\displaystyle(\ref{eq:disjunctive state chance
constraints})\xLeftarrow{(\ref{eq:boole's inequality})}$
$\displaystyle\sum^{l_{x}}_{i=1}\text{Pr}(\boldsymbol{H}_{i}\boldsymbol{x}>h_{i})\leq\epsilon_{x_{i}}$
$\displaystyle\equiv$
$\displaystyle\sum^{l_{x}}_{i=1}\text{Pr}(\boldsymbol{H}_{i}\boldsymbol{x}\leq
h_{i})\geq 1-\epsilon_{x_{i}}.$ (9)

Similarly, control joint chance constraint (6d) can be reformulated as a set
of individual chance constraints following the same arguments as before

$\sum^{l_{u}}_{i=1}\text{Pr}(\boldsymbol{G}_{i}\boldsymbol{u}\leq g_{i})\geq
1-\epsilon_{u_{i}},$ (10)

where $\epsilon_{u_{i}}=(1-\alpha_{u})/l_{u}$ is the equally distributed
control constraint risk.

######  Remark 3.1.

Allocating risk of constraint violations equally can be quite conservative
since one would preferably allocate more risk to active constraints over
inactive ones. Another approach can be optimizing for allowable violation for
each constraint as in [24], which involves higher computational complexity.

###  3.2 Deterministic Reformulation of Individual Chance Constraints

Solving the chance constraints (3.1)-(10), requires propagating the
uncertainty through the nonlinear dynamics. We adopt a linearization-based
covariance propagation as in [25][26]. Using a state-feedback control policy
$\boldsymbol{u}_{k}=\boldsymbol{v}_{k}+\boldsymbol{K}_{k}(\boldsymbol{x}_{k}-\boldsymbol{s}_{k})$,
where $\boldsymbol{K}_{k}$ are pre-stabilizing feedback gains, then the mean
and covariance of the dynamics evolve as

$\displaystyle\boldsymbol{s}_{k+1}$
$\displaystyle\approx\bar{\boldsymbol{f}}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0})+\boldsymbol{A}_{k}(\boldsymbol{s}_{k}-\boldsymbol{s}^{j}_{k})+\boldsymbol{B}_{k}(\boldsymbol{v}_{k}-\boldsymbol{v}^{j}_{k}),$
(11a) $\displaystyle\boldsymbol{\Sigma}_{\boldsymbol{x}_{k+1}}$
$\displaystyle=\boldsymbol{A}_{\text{cl}}\boldsymbol{\Sigma}_{\boldsymbol{x}_{k}}\boldsymbol{A}_{\text{cl}}^{\top}+\boldsymbol{C}_{k}\boldsymbol{\Sigma}_{\boldsymbol{\theta}}\boldsymbol{C}^{\top}_{k}+\boldsymbol{\Sigma_{w}},$
(11b)

where $\bar{\boldsymbol{f}}$ is the nominal nonlinear dynamics estimated at
current mean of the state $\boldsymbol{s}^{j}_{k}$ and controls
$\boldsymbol{v}^{j}_{k}$ of the $jth$ trajectory.
$\boldsymbol{\Sigma}_{\boldsymbol{x}_{0}}=\boldsymbol{0}$, and
$\boldsymbol{A}_{\text{cl}}\triangleq\boldsymbol{A}_{k}+\boldsymbol{B}_{k}\boldsymbol{K}_{k}$
is the closed-loop dynamics.
$\boldsymbol{A}_{k}\triangleq\frac{\partial}{\partial\boldsymbol{s}}\boldsymbol{f}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0})|_{(\boldsymbol{s}^{j}_{k},\boldsymbol{v}^{j}_{k})}$
is the Jacobian of the dynamics w.r.t. the state.
$\boldsymbol{B}_{k}\triangleq\frac{\partial}{\partial\boldsymbol{v}}\boldsymbol{f}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0})|_{(\boldsymbol{s}^{j}_{k},\boldsymbol{v}^{j}_{k})}$
is the Jacobian of the dynamics w.r.t. controls. Finally,
$\boldsymbol{C}_{k}\triangleq\frac{\partial}{\partial\boldsymbol{p}_{e,k}}\boldsymbol{f}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0})|_{(\boldsymbol{s}^{j}_{k},\boldsymbol{v}^{j}_{k})}$
represents the Jacobian of the dynamics w.r.t. the contact positions.

######  Remark 3.2.

Other approaches can be used for uncertainty propagation through nonlinear
dynamics like unscented-based transforms [27], or Generalized Polynommial
Chaos (gPC) [28]. These methods can lead to more accurate estimate of the
propagated uncertainty at the cost of significant increase in complexity.
Since computational efficiency is more important in our case (especially for
online re-planning of the trajectories), we prefer to not use these methods.

Based on Assumption (1) and the covariance propagation in (11b), we seek the
least conservative upper bounds on the state and controls individual chance-
constraints (6c)-(6d). Using the inverse of the Cumulative Density Function
(CDF) $\phi^{-1}$ of a Gaussian distribution, we arrive to a deterministic
reformulation of the chance constraints:

$\displaystyle\boldsymbol{H}_{i}\boldsymbol{s}_{k}\leq h_{i}-\eta_{x_{i,k}},$
(12a) $\displaystyle\boldsymbol{G}_{i}\boldsymbol{v}_{k}\leq
g_{i}-\eta_{u_{i,k}},$ (12b)

where
$\eta_{x_{i,k}}=\phi^{-1}(1-\epsilon_{x_{i}})\mathinner{\\!\left\lVert\boldsymbol{H}_{i}\right\rVert}_{\boldsymbol{\Sigma}_{k}}$
and
$\eta_{u_{i,k}}=\phi^{-1}(1-\epsilon_{u_{i}})\mathinner{\\!\left\lVert\boldsymbol{G}_{i}\boldsymbol{K}_{k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}$
are known as the state and control back-off bounds ensuring the satisfaction
of the individual chance constraints (6c)-(6d), respectively.

###  3.3 Deterministic Reformulation of SOCP

Given the previous reformulation of the individual chance constraints, we can
write down the following NOCP.

######  Problem 3.

NOCP with reformulated individual chance-constraints:

$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{s},\boldsymbol{v}\end{subarray}}.\,\,\,l_{f}(\boldsymbol{s}_{N})+\sum^{N-1}_{i=0}l(\boldsymbol{s}_{k},\boldsymbol{v}_{k})$
(13a) s.t.
$\displaystyle\quad\boldsymbol{s}_{k+1}=\boldsymbol{f}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0}),$
(13b) $\displaystyle\quad\boldsymbol{H}_{i,k}\boldsymbol{s}_{k}\leq
h_{i,k}-\eta_{x_{i,k}}\quad\quad\forall i\in\\{1,2,..,l_{x}\\},$ (13c)
$\displaystyle\quad\boldsymbol{G}_{i,k}\boldsymbol{v}_{k}\leq
g_{i,k}-\eta_{u_{i,k}}\quad\quad\,\forall i\in\\{1,2,..,l_{u}\\},$ (13d)
$\displaystyle\quad\boldsymbol{s}_{0}=\boldsymbol{s}(0),$ (13e)
$\displaystyle\quad\boldsymbol{s}_{f}=\boldsymbol{s}(N),$ (13f)
$\displaystyle\quad\forall k\in\\{0,1,..,N-1\\}.$ (13g)

where (13b) is now the mean of the nonlinear dynamics. In order to solve the
the above nonlinear OCP, we resort to Sequential Convex Programming (SCP),
which we explain in the next subsection.

###  3.4 SCP with L1 Trust Region Penalty Cost

SCP attempts to solve nonlinear OCPs by successively linearizing the dynamics,
costs and constraints to solve a convex sub-problem at every iteration. The
dynamics are linearized with a first-order Taylor expansion around the
previous state and control trajectories computed at the $j$-th succession.
Successive linearization introduces two well-known problems [24].

1) Artificial infeasibility: the problem becomes infeasible even if the
original nonlinear problem is feasible. The most evident example of this
arises when the problem is linearized about an unrealistically short time
horizon, so that there is no feasible control input that can satisfy the
prescribed dynamics and constraints. 2) Artificial unboundedness: the solution
takes steps far away from the validity of the linear model. In order to
mitigate artificial unboundedness, a trust-region constraint is employed.
Different approaches are adapted to tackle artificial infeasiblity. In [24],
the authors employ hard constraints and virtual controls as slack variables on
the constraints. However, [29] enforced hard constraints on the dynamics and
convex soft penalties on the rest of the constraints along with trust region
constraints. In this work we follow the same rationale as [29, 30], where the
trust region constraints $c_{i}(\boldsymbol{x})\leq 0$ are enforced as $l_{1}$
penalty cost in the form of

$\displaystyle\operatorname*{arg\,min}.\\{\operatorname*{arg\,max}.\gamma(c_{i}(\boldsymbol{x}),0)\\},$
(14a) $\displaystyle
c_{i}(\boldsymbol{x})=|\boldsymbol{x}_{k}-\boldsymbol{x}^{j}_{k}|-\Omega.$
(14b)

where $\Omega$ is the trust region radius. Notice that the above $l_{1}$
penalty cost is exact—meaning that as the penalty weight $\gamma$ gets
infinitely large, the constraint violations are driven to zero. Even though
(14a) is non-differentiable, yet it can be solved efficiently by introducing a
slack variable $t$ as follows:

$\displaystyle\min_{t}.\quad\gamma t$ (15a) s.t.
$\displaystyle\quad|\boldsymbol{x}-\boldsymbol{x}^{j}|-\Omega\leq t,$ (15b)
$\displaystyle\quad-t\leq 0.$ (15c)

In order to solve problem (3), we solve a sequence of Quadratic Programs (QPs)
in problem (4), accompanied by a trust region update mechanism based on the
accuracy ratio of the linearized model w.r.t. the nonlinear model as in
[29][26].

######  Problem 4.

Convexified QP at the $j$-th SCP iteration:

$\displaystyle\min_{\begin{subarray}{c}\boldsymbol{s},\boldsymbol{v},\boldsymbol{t}\end{subarray}}.\,\,\,l_{f}(\boldsymbol{s}_{N})+\sum^{N-1}_{i=0}l(\boldsymbol{s}_{k},\boldsymbol{v}_{k})+\gamma^{j}\sum^{N}_{i=0}t_{k}$
(16a) s.t.
$\displaystyle\quad\boldsymbol{s}_{k+1}=\bar{\boldsymbol{f}}(\boldsymbol{s}_{k},\boldsymbol{v}_{k},\boldsymbol{p}_{e,k},\boldsymbol{0})+\boldsymbol{A}_{k}(\boldsymbol{s}_{k}-\boldsymbol{s}^{j}_{k})+\boldsymbol{B}_{k}(\boldsymbol{v}_{k}-\boldsymbol{v}^{j}_{k}),$
(16b)
$\displaystyle\quad\boldsymbol{\Sigma}_{k+1}=\boldsymbol{A}_{\text{cl}}\boldsymbol{\Sigma}_{\boldsymbol{x}_{k}}\boldsymbol{A}_{\text{cl}}^{\top}+\boldsymbol{C}_{k}\boldsymbol{\Sigma}_{\boldsymbol{\theta}}\boldsymbol{C}^{\top}_{k}+\boldsymbol{\Sigma_{w}},$
(16c) $\displaystyle\quad\boldsymbol{\Sigma}_{0}=\boldsymbol{0}_{9\times 9},$
(16d) $\displaystyle\quad\boldsymbol{H}_{i,k}\boldsymbol{s}_{k}\leq
h_{i,k}-\phi^{-1}(1-\epsilon_{x_{i}})\Big{(}\mathinner{\\!\left\lVert\boldsymbol{H}_{i,k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}+\frac{\partial}{\partial{\boldsymbol{z}}}\mathinner{\\!\left\lVert\boldsymbol{H}_{i,k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}(\boldsymbol{z}_{k}-\boldsymbol{z}_{k}^{j})\Big{)},$
$\displaystyle\hskip 241.84842pt\forall i\in\\{1,2,..,l_{x}\\},$ (16e)
$\displaystyle\quad\boldsymbol{G}_{i,k}\boldsymbol{v}_{k}\leq
g_{i,k}-\phi^{-1}(1-\epsilon_{u_{i}})\Big{(}\mathinner{\\!\left\lVert\boldsymbol{G}_{i,k}\boldsymbol{K}_{k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}+\frac{\partial}{\partial{\boldsymbol{z}}}\mathinner{\\!\left\lVert\boldsymbol{G}_{i,k}\boldsymbol{K}_{k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}(\boldsymbol{z}_{k}-\boldsymbol{z}_{k}^{j})\Big{)},$
$\displaystyle\hskip 241.84842pt\forall i\in\\{1,2,..,l_{u}\\},$ (16f)
$\displaystyle\quad|\boldsymbol{\kappa}_{k}-\boldsymbol{\kappa}_{k}^{j}|-\Omega^{j}\leq
t_{k},\quad-t_{k}\leq 0,$ (16g)
$\displaystyle\quad\boldsymbol{s}_{0}=\boldsymbol{s}(0),$ (16h)
$\displaystyle\quad\boldsymbol{s}_{f}=\boldsymbol{s}(N),$ (16i)
$\displaystyle\quad\forall k\in\\{0,1,..,N-1\\}.$ (16j)

$\boldsymbol{z}_{k}\in\mathbb{R}^{9+3n_{c}}=(\boldsymbol{s}_{k},\boldsymbol{v}_{k})$
is the concatenated vector of states and controls at time $k$. Constraints
(16e)-(16f) are the linearized state and control chance constraints, where

$\displaystyle\frac{\partial}{\partial{\boldsymbol{z}}}\mathinner{\\!\left\lVert\boldsymbol{H}_{i,k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}=\frac{1}{2\mathinner{\\!\left\lVert\boldsymbol{H}_{i,k}\right\rVert}_{\boldsymbol{\Sigma}_{k}}}\Big{(}2\boldsymbol{H}^{\top}_{i}\boldsymbol{\Sigma}_{k}\frac{\partial}{\partial\boldsymbol{z}}\boldsymbol{H}_{i,k}+\sum^{n}_{i=0}\sum^{n}_{j=0}h_{i}h_{j}\frac{\partial}{\partial\boldsymbol{z}}\Sigma_{ij}\Big{)}.$
(17)

$\frac{\partial}{\partial\boldsymbol{z}}\boldsymbol{\Sigma}\in\mathbb{R}^{9\times
9\times(9+3n_{c})}$ represents the covariance derivative w.r.t.
$\boldsymbol{z}$. Notice that this term is more involved since it includes the
propagation of the tensor derivatives of the covariance matrix given the
current states and controls as well as the previous states and controls as
follows:

$\frac{\partial}{\partial\boldsymbol{z}}\boldsymbol{\Sigma}_{k+1}=\sum^{k-1}_{i=0}\boldsymbol{A}_{k}\frac{\partial}{\partial\boldsymbol{z}}\boldsymbol{\Sigma}_{k+1|i}\boldsymbol{A}^{\top}_{k}+\frac{\partial}{\partial\boldsymbol{z}}\boldsymbol{\Sigma}_{k+1|k}.$
(18)

We resort to the autodiff library JAX [31] for such computation. Finally, the
trust region constraints (16g) are enforced only on the angular momentum
$\boldsymbol{\kappa}_{k}$ since it’s the only nonlinear part in the centroidal
dynamics.

##  4 Simulations Results

In this section, we report simulation results for the quadruped robot Solo in
the Pybullet simulation environment [22]. We compare trajectories generated
using centroidal stochastic trajectory optimization against nominal trajectory
optimization for trotting and bounding gaits on challenging unknown cluttered
terrains. Offline, we warm start the centroidal SCP solver using centroidal
trajectories coming from the whole-body DDP solver Croccodyl [11]. Then, we
optimize whole-body trajectories to track back the optimized centroidal and
force trajectories from the SCP solver as illustrated pictorially in Fig. 1.
The cost weights for both whole-body DDP and centroidal SCP are summarized in
Table 1 and Table 2, respectively. Both DDP and SCP solvers were discretized
with a sampling time of $\Delta_{k}=10$ ms for a planning horizon length of
$N=165$, and motion plans were designed on a flat ground with a floor static
coefficient of friction $\mu=0.5$ for both solvers.

Figure 2: Trot motion in an unplanned cluttered environment using stochastic centroidal SCP and whole-body DDP. Table 1: Whole-body DDP cost weights. | DDP solver task weights
---|---
Task | Trot | Bound
Swing foot | 1e6 | 1e6
SCP CoM tracking | 1e3 | 1e1
SCP centroidal tracking | 1e3 | 1e3
SCP force tracking | 1e2 | 8e1
Friction cone | 2e2 | 2e0
State regulation | 1e-1 | 1e-1
Control regulation | 1e0 | 1e0
Contact impact velocity regulation | 2e1 | 2e1

During simulation (i.e. online), whole-body DDP joint-space trajectories were
tracked at a higher sampling rate of $\Delta_{k}=1$ ms using a PD control law:

$\boldsymbol{\tau}_{k}=\bar{\boldsymbol{\tau}}_{k}+\boldsymbol{K}_{p}(\boldsymbol{q}_{k}-\tilde{\boldsymbol{q}}_{k})+\boldsymbol{K}_{d}(\dot{\boldsymbol{q}}_{k}-\dot{\tilde{\boldsymbol{q}}}_{k}),$
(19)

where $\bar{\boldsymbol{\tau}}_{k}$ are the DDP optimal feedforward joint
torque controls, $\boldsymbol{q}_{k}$ and $\dot{\boldsymbol{q}}_{k}$ represent
the DDP optimal joint positions and velocities respectively. Although in
theory the optimal DDP gains could be used, it was not transferable in our
case for highly dynamic motions, especially with long horizon as the motion
diverged quickly. For that reason, we used hand-tuned PD gains for the
scenarios described in the following subsection. The chance-constraints hyper
parameters of the stochastic SCP were tuned as follows for the trotting and
bounding motions: the probability level of friction pyramid constraint
violations for every leg is $\alpha_{u}=0.1$. The covariance of the contact
position parametric uncertainties for each foot is set to
$\boldsymbol{\Sigma}_{\theta}=\text{diag}\begin{bmatrix}0.4^{2},0.4^{2},0.4^{2}\end{bmatrix}$.
The covariance of additive centroidal uncertainties is set to
$\boldsymbol{\Sigma}_{w}=\text{diag}\,\begin{bmatrix}0.85^{2},0.4^{2},0.01^{2},0.75^{2},0.4^{2},0.01^{2},0.85^{2},0.4^{2},0.01^{2}\end{bmatrix}$,
and
$\boldsymbol{\Sigma}_{w}=\text{diag}\,\begin{bmatrix}0.75^{2},0.4^{2},0.01^{2},0.85^{2},0.4^{2},0.01^{2},0.75^{2},0.4^{2},0.01^{2}\end{bmatrix}$
for the trotting and bounding motions, respectively.

###  4.1 Simulations setup

We ran a set of Monte-carlo simulations for two scenarios per motion: Scenario
1) without debris: 100 simulations on flat ground with a reduced floor
friction $\mu=0.4$, while applying random lateral force disturbances for 200
ms at the center of the robot’s base link. For trotting motion, we set
$\boldsymbol{K}_{p}=4.0*\mathbb{I}$, $\boldsymbol{K}_{d}=0.2*\mathbb{I}$. For
the bounding motion, we set $\boldsymbol{K}_{p}=3.0*\mathbb{I}$,
$\boldsymbol{K}_{d}=0.2*\mathbb{I}$. Scenario 2) with debris: 100 simulations
with reduced floor friction $\mu=0.4$, while adding unplanned debris of $2-3$
cm height ($6.6-10\%$ of the robot’s leg length) with varying orientations of
$0-17$ degrees along x and y directions as shown in Fig. 2 and Fig. 3 for
trotting and bounding motions, respectively (please refer to the video for
more details).The joint impedances were set to
$\boldsymbol{K}_{p}=5.0*\mathbb{I}$, $\boldsymbol{K}_{d}=0.2*\mathbb{I}$ for
the trotting motion, and $\boldsymbol{K}_{p}=4.7*\mathbb{I}$,
$\boldsymbol{K}_{d}=0.2*\mathbb{I}$ for the bounding motion. Further, we apply
again random lateral force impulses for 200 ms at the center of the robot’s
base.

Figure 3: Bound motion in an unplanned cluttered environment using stochastic centroidal SCP and whole-body DDP. Table 2: Centroidal SCP cost weights. | SCP solver task weights
---|---
Task | Trot | Bound
DDP CoM tracking | 1e4 | 1e4
DDP linear momentum tracking | 1e3 | 1e3
DDP angular momentum tracking | 1e5 | 1e5
Lateral force regulation per foot (x-direction) | 1e2 | 1e2
Lateral force regulation per foot (y-direction) | 1e0 | 1e2
Vertical force regulation per foot | 1e1 | 1e1
Initial trust region weight | 1e2 | 1e2

The same force disturbances were applied to the nominal and stochastic
trajectories, and were sampled from a Gaussian distribution with zero mean and
$\sigma=15$ N ($60\%$ of the robot weight). The force impulse is applied at
the same randomly sampled time instance after the first second of the motion.
We analyze the robustness of the motions generated using stochastic SCP
against their nominal counterpart by evaluating the Normalized cumulative sum
of the contact position deviations of the robot feet when a foot is in contact
with the ground (i.e. foot slippage), which reflects the saturation of the
friction pyramid constraints. The normalized cumulative sum was computed by
subtracting the average cumulative sum of the previous samples from the
current integral quantity at each point in time. Moreover, we report the
centroidal tracking performance between the generated SCP references and the
simulated trajectories.

First, we discuss the optimized contact forces generated using nominal and
stochastic SCPs, which are later tracked using whole-body DDP. In Fig. 4, we
plot the ratio between the norm of the tangential forces and the vertical
forces for dynamic trotting and bounding motions. As expected, the forces
optimized using stochastic SCP saturate less the friction cones compared to
the ones optimized with nominal SCP, especially during single support phases
where the solution of the QP is unique. This highlights the contribution of
the control back-off magnitudes, which increase along the horizon due to the
covariance propagation along the linearized dynamics (16f).

For the trotting motion (Fig. 5), trajectories designed using stochastic SCP
(our method) achieved less feet slippage mean ($26.3\%$ and $28.9\%$ for
scenario 1 and 2, respectively) than nominal SCP and an improved centroidal
tracking performance mean ($8.41\%$ and $13.0\%$). The same analysis was
carried out for the bounding motion in a more challenging terrain (Fig. 3). As
shown in Fig. 6, stochastic SCP trajectories contributed to less feet slippage
mean ($22.8\%$ and $14.8\%$ for scenarios 1 and 2, respectively) than nominal
SCP, and an improved centroidal tracking performance mean ($25.6\%$ and
$13.6\%$).

Figure 4: Ratio of norm of tangential forces w.r.t. vertical force for a
trotting motion (left) and a bounding motion (right).

##  5 Discussion and Conclusions

In this work, we used nonlinear stochastic trajectory optimization for
generating robust centroidal momentum trajectories for legged robots that take
into account additive uncertainties on the centroidal dynamics as well as
parametric uncertainties on the contact positions. We used a linearization-
based covariance propagation for resolving the stochastic nonlinear dynamics.
Furthermore, we resolved the friction pyramid joint chance constraints by
designing proper upper bounds (back-offs) at each point in time on the
individual hyper-planes forming the friction pyramid polytopes. Finally, we
presented a whole-body trajectory optimization framework that alternates
between stochastic centroidal trajectory optimization and whole-body
trajectory optimization for generating feasible robust whole-body motions. We
used our framework to generate trotting and bounding dynamic gaits for the
quadruped robot Solo. We then tracked these trajectories in a Pybullet physics
simulator, introducing different disturbance realizations and contact
uncertainties. The results show that our approach generated safer motions by
contributing to less average contact slippage, as well as improved centroidal
tracking performance over deterministic trajectory optimization. Although the
current stochastic SCP approach does not require additional optimization
variables over a deterministic approach SCP, the computational complexity is
relatively higher due to uncertainty propagation and the additional tensor
derivatives required for solving the linearized chance-constraints.

Another limitation of the current stochastic SCP approach lies in the accuracy
of uncertainties propagation through the linearized dynamics, which might be
hindered for long horizons. However, we believe that this might not be an
issue in practice when applied in receding horizon fashion. To this end, we
plan to extend the current framework to nonlinear stochastic MPC in the
future.

Figure 5: Normalized cumulative sum of feet slippage norm (left) and
centroidal tracking cost (right) for a trotting motion.

Figure 6: Normalized cumulative sum of feet slippage norm (left) and
centroidal tracking cost (right) for a bounding motion.

## References

  * [1] Mordatch, I., Todorov, E., Popović, Z.: Discovery of complex behaviors through contact-invariant optimization. ACM Trans. Graph. 31(4) (jul 2012)
  * [2] Posa, M., Cantu, C., Tedrake, R.: A direct method for trajectory optimization of rigid bodies through contact. The International Journal of Robotics Research 33(1), 69–81 (2014)
  * [3] Winkler, A.W., Bellicoso, C.D., Hutter, M., Buchli, J.: Gait and trajectory optimization for legged systems through phase-based end-effector parameterization. IEEE Robotics and Automation Letters 3(3), 1560–1567 (2018)
  * [4] Carpentier, J., Wieber, P.B.: Recent progress in legged robots locomotion control. Current Robotics Reports 2(3), 231–238 (2021)
  * [5] Sleiman, J.P., Farshidian, F., Minniti, M.V., Hutter, M.: A unified mpc framework for whole-body dynamic locomotion and manipulation. IEEE Robotics and Automation Letters 6(3), 4688–4695 (2021)
  * [6] Meduri, A., Shah, P., Viereck, J., Khadiv, M., Havoutis, I., Righetti, L.: Biconmp: A nonlinear model predictive control framework for whole body motion planning. CoRR abs/2201.07601 (2022), https://arxiv.org/abs/2201.07601
  * [7] Villa, N.A., Wieber, P.B.: Model predictive control of biped walking with bounded uncertainties. In: 2017 IEEE-RAS 17th International Conference on Humanoid Robotics (Humanoids). pp. 836–841 (2017)
  * [8] Gazar, A., Khadiv, M., Prete, A.D., Righetti, L.: Stochastic and robust mpc for bipedal locomotion: A comparative study on robustness and performance. In: 2020 IEEE-RAS 20th International Conference on Humanoid Robots (Humanoids). pp. 61–68 (2021)
  * [9] Yeganegi, M.H., Khadiv, M., Moosavian, S.A.A., Zhu, J.J., Del Prete, A., Righetti, L.: Robust humanoid locomotion using trajectory optimization and sample-efficient learning. In: 2019 IEEE-RAS 19th International Conference on Humanoid Robots (Humanoids). pp. 170–177 (2019)
  * [10] Tassa, Y., Erez, T., Todorov, E.: Synthesis and stabilization of complex behaviors through online trajectory optimization. In: 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems. pp. 4906–4913 (2012)
  * [11] Mastalli, C., Budhiraja, R., Merkt, W., Saurel, G., Hammoud, B., Naveau, M., Carpentier, J., Righetti, L., Vijayakumar, S., Mansard, N.: Crocoddyl: An efficient and versatile framework for multi-contact optimal control. In: 2020 IEEE International Conference on Robotics and Automation (ICRA). pp. 2536–2542 (2020)
  * [12] Mordatch, I., Lowrey, K., Todorov, E.: Ensemble-cio: Full-body dynamic motion planning that transfers to physical humanoids. In: 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). pp. 5307–5314 (2015)
  * [13] Hammoud, B., Khadiv, M., Righetti, L.: Impedance optimization for uncertain contact interactions through risk sensitive optimal control. IEEE Robotics and Automation Letters 6(3), 4766–4773 (2021)
  * [14] Drnach, L., Zhang, J.Z., Zhao, Y.: Mediating between contact feasibility and robustness of trajectory optimization through chance complementarity constraints. Frontiers in Robotics and AI 8 (2021)
  * [15] Orin, D.E., Goswami, A.: Centroidal momentum matrix of a humanoid robot: Structure and properties. In: 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems. pp. 653–659 (2008)
  * [16] Wieber, P.B., Tedrake, R., Kuindersma, S.: Modeling and Control of Legged Robots, pp. 1203–1234. Springer International Publishing, Cham (2016)
  * [17] Herzog, A., Schaal, S., Righetti, L.: Structured contact force optimization for kino-dynamic motion generation. In: 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). pp. 2703–2710 (2016)
  * [18] Brasseur, C., Sherikov, A., Collette, C., Dimitrov, D., Wieber, P.B.: A robust linear mpc approach to online generation of 3d biped walking motion. In: 2015 IEEE-RAS 15th International Conference on Humanoid Robots (Humanoids). pp. 595–601 (2015)
  * [19] Dai, H., Tedrake, R.: Planning robust walking motion on uneven terrain via convex optimization. In: 2016 IEEE-RAS 16th International Conference on Humanoid Robots (Humanoids). pp. 579–586 (2016)
  * [20] Budhiraja, R., Carpentier, J., Mansard, N.: Dynamics consensus between centroidal and whole-body models for locomotion of legged robots. In: 2019 International Conference on Robotics and Automation (ICRA). pp. 6727–6733 (2019)
  * [21] Grimminger, F., Meduri, A., Khadiv, M., Viereck, J., Wüthrich, M., Naveau, M., Berenz, V., Heim, S., Widmaier, F., Flayols, T., et al.: An open torque-controlled modular robot architecture for legged locomotion research. IEEE Robotics and Automation Letters 5(2) (2020)
  * [22] Coumans, E., Bai, Y.: Pybullet, a python module for physics simulation for games, robotics and machine learning. http://pybullet.org (2016–2021)
  * [23] Ono, M., Blackmore, L., Williams, B.C.: Chance constrained finite horizon optimal control with nonconvex constraints. In: Proceedings of the 2010 American Control Conference. pp. 1145–1152 (2010)
  * [24] Ma, Y., Vichik, S., Borrelli, F.: Fast stochastic mpc with optimal risk allocation applied to building control systems. In: 2012 IEEE 51st IEEE Conference on Decision and Control (CDC). pp. 7559–7564 (2012)
  * [25] Zhu, H., Alonso-Mora, J.: Chance-constrained collision avoidance for mavs in dynamic environments. IEEE Robotics and Automation Letters 4(2), 776–783 (2019)
  * [26] Lew, T., Bonalli, R., Pavone, M.: Chance-constrained sequential convex programming for robust trajectory optimization. In: 2020 European Control Conference (ECC) (2020)
  * [27] Plancher, B., Manchester, Z., Kuindersma, S.: Constrained unscented dynamic programming. In: 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). pp. 5674–5680 (2017)
  * [28] Nakka, Y.K., Chung, S.J.: Trajectory optimization for chance-constrained nonlinear stochastic systems. In: 2019 IEEE 58th Conference on Decision and Control (CDC). pp. 3811–3818 (2019)
  * [29] Bonalli, R., Cauligi, A., Bylard, A., Pavone, M.: Gusto: Guaranteed sequential trajectory optimization via sequential convex programming. In: 2019 International Conference on Robotics and Automation (ICRA). pp. 6741–6747 (2019)
  * [30] Schulman, J., Duan, Y., Ho, J., Lee, A., Awwal, I., Bradlow, H., Pan, J., Patil, S., Goldberg, K., Abbeel, P.: Motion planning with sequential convex optimization and convex collision checking. The International Journal of Robotics Research (2014)
  * [31] Bradbury, J., Frostig, R., Hawkins, P., Johnson, M.J., Leary, C., Maclaurin, D., Necula, G., Paszke, A., VanderPlas, J., Wanderman-Milne, S., Zhang, Q.: JAX: composable transformations of Python+NumPy programs, http://github.com/google/jax


# Kählerity of Einstein four-manifolds

Xiaolong Li  Department of Mathematics, Statistics and Physics, Wichita State
University, Wichita, KS, 67260<EMAIL_ADDRESS>and  Yongjia Zhang
School of Mathematical Sciences, Shanghai Jiao Tong University, Shanghai,
China, 200240<EMAIL_ADDRESS>
###### Abstract.

We prove that a closed oriented Einstein four-manifold is either anti-self-
dual or (after passing to a double Riemannian cover if necessary) Kähler-
Einstein, provided that $\lambda_{2}\geq-\frac{S}{12}$, where $\lambda_{2}$ is
the middle eigenvalue of the self-dual Weyl tensor $W^{+}$ and $S$ is the
scalar curvature. An analogous result holds for closed oriented four-manifolds
with $\delta W^{+}=0$.

###### Key words and phrases:

Einstein four-manifolds, Kähler-Einstein, self-dual Weyl tensor, curvature
restriction

###### 2020 Mathematics Subject Classification:

53C21, 53C25, 53C55

The first author’s research is partially supported by a start-up grant at
Wichita State University

##  1\. Introduction

We are concerned with the following question in this paper.

######  Question 1.1.

Under what curvature conditions is a closed Einstein four-manifold Kähler?

There have been several answers to this question. Derdzinski proved a
fundamental result [Der83] that if the self-dual Weyl tensor $W^{+}$ of a
closed oriented Einstein four-manifold is parallel and has at most two
distinct eigenvalues at every point, then, after passing to a double
Riemannian cover if necessary, the metric is Kähler. A result of Micallef and
Wang [MW93, Theorem 4.2] asserts that a closed oriented four-manifold with
harmonic self-dual Weyl tensor ($\delta W^{+}=0$) and half nonnegative
isotropic curvature either is anti-self-dual or has a Kähler universal cover
with constant positive scalar curvature (see also [GL99], [RS16], [FKP14], and
[Wu17] for alternative proofs).

The purpose of this paper is to provide another curvature condition that
characterizes closed Kähler-Einstein four-manifolds. To state our results, let
us recall that the bundle $\Lambda^{2}$ of two-forms over an oriented (that
is, orientable with an appointed orientation) Riemannian four-manifold
$(M^{4},g)$ admits an orthogonal decomposition

$\Lambda^{2}=\Lambda^{+}\oplus\Lambda^{-},$

with $\Lambda^{+}$ and $\Lambda^{-}$ being the $+1$ and $-1$ eigenspaces of
the Hodge star operator, respectively. Sections of $\Lambda^{+}$ are called
self-dual two-forms, while sections of $\Lambda^{-}$ are called anti-self-dual
two-forms. The Riemann curvature tensor may be identified with a self-adjoint
linear map

$\operatorname{Rm}:\Lambda^{2}\to\Lambda^{2}$

called the curvature operator, via

$\operatorname{Rm}(e_{i}\wedge e_{j})=\frac{1}{2}R_{ijkl}e_{k}\wedge e_{l}.$

Corresponding to the decomposition $\Lambda^{2}=\Lambda^{+}\oplus\Lambda^{-}$,
the Riemann curvature operator $\operatorname{Rm}$ can be decomposed into the
following irreducible pieces:

$\displaystyle\operatorname{Rm}=\left[\begin{array}[]{c|c}&\\\
W^{+}+\frac{S}{12}I&\overset{\circ}{\operatorname{Ric}}\\\ &\\\ \hline\cr&\\\
\overset{\circ}{\operatorname{Ric}}&W^{-}+\frac{S}{12}I\\\
&\end{array}\right],$

where $S$ is the scalar curvature, $I$ is the identity map,
$\overset{\circ}{\operatorname{Ric}}=\operatorname{Ric}-\tfrac{S}{4}g$ is the
traceless Ricci tensor, and $W^{\pm}:\Lambda^{\pm}\to\Lambda^{\pm}$ are called
the self-dual and anti-self-dual parts of the Weyl curvature operator
$W:\Lambda^{2}\to\Lambda^{2}$. It is also convenient to denote by
$\operatorname{Rm}^{+}$ and $\operatorname{Rm}^{-}$ the restriction of the
Riemann curvature operator on $\Lambda^{+}$ and $\Lambda^{-}$, respectively.
In other words, $\operatorname{Rm}^{\pm}=W^{\pm}+\frac{S}{12}I$. A four-
manifold $(M^{4},g)$ is said to be anti-self-dual if $W^{+}\equiv 0$ and self-
dual if $W^{-}\equiv 0$.

It is well-known that on any Kähler surface $(M^{4},J,g)$ with the natural
orientation (in the sense that the Kähler form is self-dual), the self-dual
Weyl operator $W^{+}:\Lambda^{+}\to\Lambda^{+}$ is given by

$\displaystyle W^{+}=\left[\begin{array}[]{ccc}\tfrac{S}{6}&&\\\
&-\tfrac{S}{12}&\\\ &&-\tfrac{S}{12}\end{array}\right].$

In particular, if $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ are the
eigenvalues of $W^{+}$, then we have

$\lambda_{2}=-\frac{S}{12}$

on any Kähler surface, regardless of the sign of $S$. Equivalently, the middle
eigenvalue of $\operatorname{Rm}^{+}$ vanishes on any Kähler surfaces.

Our main result states that if a closed oriented Einstein four-manifold
satisfies the condition

$\lambda_{2}\geq-\frac{S}{12},$

then it is either anti-self-dual or Kähler-Einstein (after passing to a double
cover if necessary). More precisely, we prove

######  Theorem 1.1.

Let $(M^{4},g)$ be a closed oriented Einstein four-manifold. Denote by
$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ the eigenvalues of $W^{+}$ and by
$S$ the scalar curvature. If $\lambda_{2}\geq-\frac{S}{12}$ everywhere, then
one of the following assertions holds:

  1. (1)

$S\neq 0$ and $(M,g)$ is Kähler-Einstein, or has a double Riemannian cover
$\pi:(M^{*},g^{*})\to(M,g)$, where $g^{*}=\pi^{*}g$, such that $(M^{*},g^{*})$
is Kähler-Einstein.

  2. (2)

$(M,g)$ is anti-self-dual with nonnegative scalar curvature. Moreover,

     * •

if $S>0$, then $M$ is isometric to a quotient of $\mathbb{S}^{4}$ with the
round metric or $\mathbb{CP}^{2}$ with the Fubini-Study metric.

     * •

if $S=0$, then either $M$ is flat or its universal cover is a K3 surface with
the Calabi-Yau metric.

In the nonnegative scalar curvature case, Theorem 1.1 provides a
generalization of the result of Micallef and Wang [MW93], as the condition
$\lambda_{2}\geq-\frac{S}{12}$ is weaker than half nonnegative isotropic
curvature. Recall that an oriented four-manifold is said to have half
nonnegative (resp. positive) isotropic curvature if $\operatorname{Rm}^{+}$ is
two-nonnegative (resp. two-positive), i.e., the sum of the smallest two
eigenvalues of $\operatorname{Rm}^{+}$ is nonnegative (resp. positive). In
view of $\operatorname{Rm}^{+}=W^{+}+\frac{S}{12}I$, we see that half
nonnegative isotropic curvature is equivalent to

$S\geq 0\quad\text{ and }\quad\lambda_{1}+\lambda_{2}\geq-\tfrac{S}{6},$

so it is stronger than the condition $\lambda_{2}\geq-\frac{S}{12}$.

In the negative scalar curvature case, Theorem 1.1, to the best of our
knowledge, seems to be the only answer to Question 1.1 other than that of
Derdzinski [Der83].

As a corollary of Theorem 1.1, we have a sufficient and necessary condition
for a closed Einstein four-manifold to be Kähler.

######  Corollary 1.2.

Let $(M^{4},g)$ be a closed and oriented (positive or negative) Einstein four-
manifold which satisfies $H^{1}(M;\mathbb{Z}_{2})=0$ and is not isometric to
$\mathbb{S}^{4}$ with the round metric. Then $(M^{4},g)$ is Kähler if and only
if $\lambda_{2}(W^{+})\geq-\frac{S}{12}$ everywhere.

In addition, Theorem 1.1 also implies some interesting classification results.
In the positive scalar curvature case, we have

######  Corollary 1.3.

Let $(M^{4},g)$ be a closed simply-connected Einstein four-manifold with
positive scalar curvature. If $\lambda_{2}(W^{+})>-\frac{S}{12}$ on $M$, then
$W^{+}\equiv 0$ and $M$ is isometric to $\mathbb{S}^{4}$ with the round metric
or $\mathbb{CP}^{2}$ with the Fubini-Study metric.

######  Remark 1.4.

Corollary 1.3 was proved by Polombo [Pol92] assuming the stronger condition of
half positive isotropic curvature, which is equivalent to

$S>0\quad\text{ and }\quad\lambda_{1}+\lambda_{2}>-\tfrac{S}{6}.$

Alternative proofs of Polombo’s result can be found in [MW93], [FKP14],
[RS16], and [Wu17].

######  Remark 1.5.

The assumption $\lambda_{2}>-\frac{S}{12}$ in Corollary 1.3 cannot be relaxed
in view of the Einstein manifold $\mathbb{S}^{2}\times\mathbb{S}^{2}$, which
has $\lambda_{2}=-\frac{S}{12}$ everywhere.

In the Ricci-flat case, we have

######  Corollary 1.6.

Let $(M^{4},g)$ be a closed oriented Ricci-flat four-manifold. If
$\lambda_{2}(W^{+})\geq 0$ on $M$, then $W^{+}\equiv 0$, and either $M$ is
flat or its universal cover is a K3 surface with the Calabi-Yau metric.

######  Remark 1.7.

Corollary 1.6 was proved by Micallef and Wang [MW93] under the stronger
condition of half-nonnegative isotropic curvature. Alternative proofs of their
result can be found in [FKP14], [RS16], and [Wu17]. However, it seems that
none of these proofs work under our weaker curvature assumption.

With a bit more efforts, the method employed to prove Theorem 1.1 also yields
analogous results for closed four-manifolds with harmonic (namely, divergence
free) self-dual Weyl tensor, i.e.

(1.1) $\delta W^{+}:=-\nabla\cdot W^{+}=0.$

All Einstein metrics satisfy (1.1) as a consequence of the second Bianchi
identity. However, (1.1) is actually much weaker than the Einstein condition.

######  Theorem 1.8.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$.
Denote by $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ the eigenvalues of
$W^{+}$ and by $S$ the scalar curvature. If $\lambda_{2}\geq-\frac{S}{12}$
everywhere, then one of the following assertions holds:

  1. (1)

$(M^{4},g)$ has constant nonzero scalar curvature, and $(M^{4},g)$ is Kähler,
or has a double Riemannian cover $\pi:(M^{*},g^{*})\to(M,g)$, where
$g^{*}=\pi^{*}g$, such that $(M^{*},g^{*})$ is Kähler.

  2. (2)

$(M^{4},g)$ is anti-self-dual with nonnegative scalar curvature.

We would like to point out that both Theorem 1.1 and Theorem 1.8 are also
motivated by an important work of Peng Wu [Wu21], in which he found a
beautiful characterization of conformally Kähler Einstein metrics of positive
scalar curvature on closed oriented four-manifolds via the condition
$\det(W^{+})>0$. LeBrun [LeB21], based on his earlier work in [LeB15],
provided an entirely different proof of the result and relaxed the assumptions
to $W^{+}\neq 0$ and
$|W^{+}|^{-3}\det(W^{+})\geq-\tfrac{5\sqrt{2}}{21\sqrt{21}}$. In comparison
with Wu and Lebrun’s result, Theorem 1.1 has the advantage of giving a
curvature characterization in the negative scalar curvature case, but has the
disadvantage of using a closed curvature condition. Moreover, the condition
$\lambda_{2}\geq-\frac{S}{12}$ is not conformally invariant, thus it does not
allow any conformally Kähler manifolds.

Theorem 1.8 implies the following results, which are known under the stronger
assumption of half positive (resp. nonnegative) isotropic curvature (see
[Pol92], [MW93], [FKP14], [RS16], and [Wu17]).

######  Corollary 1.9.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$ and
positive scalar curvature. If $\lambda_{2}>-\frac{S}{12}$ on $M$, then
$W^{+}\equiv 0$.

######  Corollary 1.10.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$ and
nonnonegative scalar curvature. If $\lambda_{2}\geq 0$ on $M$, then
$W^{+}\equiv 0$.

To conclude this introduction, we give a brief discussion of our strategies to
prove the above-mentioned results. In order to prove Theorem 1.1, we apply the
maximum principle to the following partial differential inequality

(1.2) $\Delta(\lambda_{3}-\lambda_{1})\geq
6(\lambda_{3}-\lambda_{1})(\lambda_{2}+\tfrac{S}{12})$

to conclude that $\lambda_{3}-\lambda_{1}$ must be a constant. Moreover, we
must have $\lambda_{2}=-\frac{S}{12}$ everywhere unless $M$ is anti-self-dual.
This together with the Einstein condition implies that all the $\lambda_{i}$’s
are constant functions, and their values can be read from the differential
inequalities satisfied by them. The desired Kählerity then follows from the
work of Derdzinski [Der83].

The key to prove Theorem 1.8 is to show that the scalar curvature must be
constant unless $M$ is anti-self-dual. We achieve this by picking up some
extra gradient terms in (1.2), which are used to conclude that the
$\lambda_{i}$’s are constants. It is worth mentioning that these helpful extra
terms are used in Wu [Wu21] as well.

The proofs of Theorem 1.1 and Theorem 1.8 are given in Section 2 and Section
3, respectively.

##  2\. The Einstein Case

In this section, we prove Theorem 1.1. First of all, we prove two technical
propositions which will also be applied in Section 3. Note that Proposition
2.1 assumes only harmonic self-dual Weyl tensor, while Proposition 2.2
assumes, in addition, constant scalar curvature; these assumptions are
consequences of the Einstein condition.

######  Proposition 2.1.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$.
Denote by $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ the eigenvalues of
$W^{+}$. If $\lambda_{2}\geq-\frac{S}{12}$ everywhere, then either
$W^{+}\equiv 0$ or $\lambda_{3}-\lambda_{1}$ is equal to a positive constant
and $\lambda_{2}=-\frac{S}{12}$ everywhere.

######  Proposition 2.2.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$ and
with constant scalar curvature $S$. Denote by
$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ the eigenvalues of $W^{+}$. If
$\lambda_{2}\geq-\frac{S}{12}$ everywhere and $W^{+}\not\equiv 0$, then the
following statements hold:

  1. (1)

$\lambda_{2}=-\frac{S}{12}$;

  2. (2)

$\lambda_{1}=-\frac{S}{12}$ and $\lambda_{3}=\frac{S}{6}$ if $S>0$, and
$\lambda_{1}=\frac{S}{6}$ and $\lambda_{3}=-\frac{S}{12}$ if $S<0$;

  3. (3)

$\nabla W^{+}=0$.

###### Proof of Proposition 2.1.

According to [MW93, page 664], $\delta W^{+}=0$ is equivalent to the
Weitzenböck formula

(2.1) $\Delta W^{+}=\tfrac{S}{2}W^{+}-2(W^{+})^{2}-4(W^{+})^{\\#},$

where the $(W^{+})^{\\#}$ is the adjoint matrix of $W^{+}$. The reader may
consult [Ham86] or [CLN06] for more information. If we choose an orthonormal
basis $\\{\omega_{i}\\}_{i=1}^{3}$ of $\Lambda^{+}$ that diagonalizes $W^{+}$
with eigenvalues $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$, then, with
respect to this basis, we have

$\displaystyle(W^{+})^{2}=\left[\begin{array}[]{ccc}\lambda_{1}^{2}&&\\\
&\lambda_{2}^{2}&\\\ &&\lambda_{3}^{2}\end{array}\right]$

and

$\displaystyle(W^{+})^{\\#}=\left[\begin{array}[]{ccc}\lambda_{2}\lambda_{3}&&\\\
&\lambda_{1}\lambda_{2}&\\\ &&\lambda_{1}\lambda_{2}\end{array}\right].$

Therefore, the Weitzenböck formula (2.1) implies the following differential
inequalities for the eigenvalues $\lambda_{1}$ and $\lambda_{3}$

(2.2) $\displaystyle\Delta\lambda_{1}$ $\displaystyle\leq$
$\displaystyle\tfrac{S}{2}\lambda_{1}-2\lambda_{1}^{2}-4\lambda_{2}\lambda_{3},$
(2.3) $\displaystyle\Delta\lambda_{3}$ $\displaystyle\geq$
$\displaystyle\tfrac{S}{2}\lambda_{3}-2\lambda_{3}^{2}-4\lambda_{1}\lambda_{2}.$

Since $\lambda_{1}$ and $\lambda_{3}$ are only locally Lipschitz functions on
$M$, the inequalities above are all understood in the sense of barrier (see
[Cal58]). To prove (2.2) and (2.3), one may construct the barriers as follows.
Fix $p\in M$ and let $\omega\in\Lambda^{+}_{p}$ be the unit eigenvector of
$W^{+}_{p}$ associated with $\lambda_{3}$. Now we may extend $\omega$ to a
neighborhood of $p$ by parallel transport and compute using the normal
coordinates at $p$:

$\displaystyle\Delta(W^{+}(\omega,\omega))$ $\displaystyle=(\Delta
W^{+})(\omega,\omega)$
$\displaystyle=\tfrac{S}{2}W^{+}(\omega,\omega)-2|W^{+}(\omega)|^{2}-4(W^{+})^{\\#}(\omega,\omega)$
$\displaystyle=\tfrac{S}{2}\lambda_{3}-2\lambda_{3}^{2}-4\lambda_{1}\lambda_{2}\quad\text{
at }p.$

Obviously, $W^{+}(\omega,\omega)$ is a lower barrier of $\lambda_{3}$ at $p$,
and this proves (2.3); the proof of (2.2) is similar.

Subtracting (2.2) from (2.3) yields

$\displaystyle\Delta(\lambda_{3}-\lambda_{1})$ $\displaystyle\geq$
$\displaystyle\tfrac{S}{2}(\lambda_{3}-\lambda_{1})-2(\lambda_{3}^{2}-\lambda_{1}^{2})+4\lambda_{2}(\lambda_{3}-\lambda_{1})$
$\displaystyle=$
$\displaystyle(\lambda_{3}-\lambda_{1})\left(\tfrac{S}{2}-2\lambda_{1}-2\lambda_{3}+4\lambda_{2}\right)$
$\displaystyle=$ $\displaystyle
6(\lambda_{3}-\lambda_{1})\left(\lambda_{2}+\tfrac{S}{12}\right),$

where we have used $\lambda_{1}+\lambda_{2}+\lambda_{3}=0$ in the last step.
It follows from $\lambda_{2}\geq-\frac{S}{12}$ that

$\Delta(\lambda_{3}-\lambda_{1})\geq
6(\lambda_{3}-\lambda_{1})\left(\lambda_{2}+\tfrac{S}{12}\right)\geq 0,$

in the sense of barrier. By the maximum principle, we conclude that
$\lambda_{3}-\lambda_{1}\equiv c\geq 0$. Note that $c=0$ implies $W^{+}\equiv
0$. In case $c>0$, the desired equality $\lambda_{2}=-\frac{S}{12}$ follows
from

$0=\Delta(\lambda_{3}-\lambda_{1})\geq
6(\lambda_{3}-\lambda_{1})\left(\lambda_{2}+\tfrac{S}{12}\right).$

The proof is finished. $\square$

###### Proof of Proposition 2.2.

Part (1) follows immediately from Proposition 2.1 since $W^{+}\not\equiv 0$.

To prove part (2), we first notice that $\lambda_{3}-\lambda_{1}$ is equal to
a positive constant by Proposition 2.1. Now the constant scalar curvature
assumption implies that $\lambda_{1}+\lambda_{3}$ is a constant function, as

$\lambda_{1}+\lambda_{3}=-\lambda_{2}=\tfrac{S}{12}.$

It then follows that both $\lambda_{1}$ and $\lambda_{3}$ must be constant
functions. Substituting $\lambda_{2}=-\frac{S}{12}$ and
$\lambda_{1}+\lambda_{2}+\lambda_{3}=0$ into the differential inequalities
satisfied by $\lambda_{1}$ and $\lambda_{3}$, we obtain that

$\displaystyle 0=\Delta\lambda_{1}$ $\displaystyle\leq$
$\displaystyle\tfrac{S}{2}\lambda_{1}-2\lambda_{1}^{2}-4\lambda_{2}\lambda_{3}=-2(\lambda_{1}+\tfrac{S}{12})(\lambda_{1}-\tfrac{S}{6})$
$\displaystyle 0=\Delta\lambda_{3}$ $\displaystyle\geq$
$\displaystyle\tfrac{S}{2}\lambda_{3}-2\lambda_{3}^{2}-4\lambda_{1}\lambda_{2}=-2(\lambda_{3}+\tfrac{S}{12})(\lambda_{3}-\tfrac{S}{6}).$

One easily reads from above inequalities that we must have

$\displaystyle\lambda_{1}=\lambda_{3}=0,$ $\displaystyle\text{ if }S=0;$
$\displaystyle\lambda_{1}=-\tfrac{S}{12}\text{ and }\lambda_{3}=\tfrac{S}{6},$
$\displaystyle\text{ if }S>0;$ $\displaystyle\lambda_{1}=\tfrac{S}{6}\text{
and }\lambda_{3}=-\tfrac{S}{12},$ $\displaystyle\text{ if }S<0.$

Now part (2) is proved.

Part (3) follows immediately from part (2) and the constant scalar curvature
assumption.

$\square$

We now give the proof of Theorem 1.1.

###### Proof of Theorem 1.1.

If $M$ is anti-self-dual, then $0=\lambda_{2}\geq-S/12$ implies that $S$ is a
nonnegative constant. By Hitchin’s classification [Hit74] of half conformally
flat Einstein four-manifolds with nonnegative scalar curvature (see also
[Bes08, Theorem 13.30]), $(M^{4},g)$ is one of the manifolds as described in
part (2) of the statement.

If $M$ is not anti-self-dual, then by Proposition 2.2, we have that $\nabla
W^{+}=0$ and $W^{+}$ has at most two distinct eigenvalues at every point. One
can then invoke the result of Derdzinski [Der83, Therorem 2] to conclude that
either $M$ is Kähler or it has a double cover that is Kähler. Below we provide
a more direct proof using an elegant argument of Lebrun [LeB21], which is also
applicable to the proof of Theorem 1.8. We shall only present the case where
$M$ has negative constant scalar curvature, as it can be easily adapted to the
positive scalar curvature case by flipping signs and reversing directions of
the related inequalities.

Since $\lambda_{1}=S/6$ is an isolated eigenvalue of $W^{+}$, we have that the
corresponding eigenspaces of $W^{+}$ at each point on $M$ form a one-
dimensional subbundle of $\Lambda^{+}$. Denote this bundle by $L$. If $L$ is a
trivial bundle (in particular, if $H^{1}(M;\mathbb{Z}_{2})=0$), then we can
find a non-vanishing global section $\omega$ of $L$ with constant norm
everywhere. If $L$ does not have a non-vanishing global section, then we may
let $M^{*}$ be all the elements in $L$ with constant unit norm, and the
restriction of the bundle projection $\pi:M^{*}\to M$ is a double cover. If we
equip $M^{*}$ with the Remannian metric $g^{*}=\pi^{*}g$, then $(M^{*},g^{*})$
also satisfies Proposition 2.1 and $W^{+}$ admits a global eigenvector section
associated with $\lambda_{1}=S/6$. Henceforth, we will assume that $L$ admits
a global section $\omega$ with constant norm $|\omega|\equiv\sqrt{2}$.

Since $\omega$, when viewed as an operator on $TM$, satisfies

$\omega\cdot\omega=-\operatorname{id},$

we have that $\omega$ provides an almost complex structure. We need only to
show that $\omega$ is a closed form.

We first observe that

(2.4) $W^{+}(\nabla_{X}\omega,\nabla_{X}\omega)\geq 0.$

for any tangent vector $X$. This is because $\omega$ has constant norm, which
implies that $\langle\nabla_{X}\omega,\omega\rangle=0$ and hence
$\nabla_{X}\omega$ is in the eigenspace of $\lambda_{2}=\lambda_{3}=-S/12>0$.

Now recall that when $\delta W^{+}=0$, the following Weitzenböck formula holds
(c.f. [MW93, (4.1c)] or [LeB21, Equation (9)])

(2.5) $\displaystyle\Delta W^{+}=\tfrac{S}{2}W^{+}-6W^{+}\circ
W^{+}+2|W^{+}|^{2}I.$

It is easy to check that (2.5) is equivalent to (2.1). We now take inner
product of (2.5) and $\omega\otimes\omega$, and integrate over $M$ to get

$\displaystyle 0$ $\displaystyle=\int_{M}\left\langle-\Delta
W^{+}+\tfrac{S}{2}W^{+}-6W^{+}\circ
W^{+}+2|W^{+}|^{2}I,\omega\otimes\omega\right\rangle dg$
$\displaystyle=\int_{M}\Big{(}-2W^{+}(\omega,\Delta\omega)-2W^{+}(\nabla_{e}\omega,\nabla^{e}\omega)$
$\displaystyle\qquad+\tfrac{S}{2}\lambda_{1}|\omega|^{2}-6\lambda_{1}^{2}|\omega|^{2}+2(\lambda_{1}^{2}+\lambda_{2}^{2}+\lambda_{3}^{2})|\omega|^{2}\Big{)}dg$
$\displaystyle\leq-2\int_{M}\lambda_{1}\langle\omega,\Delta\omega\rangle dg$
$\displaystyle=\tfrac{S}{3}\int_{M}|\nabla\omega|^{2}dg,$

where we have applied Proposition 2.2(1)(2) and (2.4). Since $S$ is a negative
constant, we have that $\nabla\omega$ vanishes everywhere. It follows that
$\omega$ is a Kähler form.

In case where $W^{+}\not\equiv 0$ and $S$ is a positive constant, we can also
find a Kähler form in the same way.

$\square$

It is clear that Corollary 1.2, Corollary 1.3, and Corollary 1.6 can be easily
observed from the proof of Theorem 1.1.

##  3\. The Harmonic self-dual Weyl case

If $M$ is only assumed to have harmonic self-dual Weyl tensor, then the proof
given in the previous section breaks down because we cannot conclude that both
$\lambda_{1}$ and $\lambda_{3}$ are constant functions without the constant
scalar curvature assumption. We shall show, by improving the partial
differential inequality satisfied by $\lambda_{3}-\lambda_{1}$, that the
scalar curvature must be a constant unless $W^{+}\equiv 0$, as a consequence
of $\delta W^{+}=0$ and $\lambda_{2}\geq-\frac{S}{12}$.

######  Proposition 3.1.

Let $(M^{4},g)$ be a closed oriented four-manifold with $\delta W^{+}=0$.
Denote by $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ the eigenvalues of
$W^{+}$. If $\lambda_{2}\geq-\frac{S}{12}$ everywhere, then either $M$ is
anti-self-dual with nonnegative scalar curvature, or $M$ has constant scalar
curvature.

###### Proof.

Let us recall some computations in [Der83]. For any $x\in M$, we can choose an
orthogonal basis $\omega_{1}$, $\omega_{2}$, $\omega_{3}$ of
$\Lambda^{+}_{x}$, consisting of eigenvectors of $W^{+}$ such that

(3.1) $|\omega_{1}|^{2}=|\omega_{2}|^{2}=|\omega_{3}|^{2}=2.$

Consequently, we have that, at $x$

(3.2)
$W^{+}=\frac{1}{2}\left(\lambda_{1}\omega_{1}\otimes\omega_{1}+\lambda_{2}\omega_{2}\otimes\omega_{2}+\lambda_{3}\omega_{3}\otimes\omega_{3}\right)$

with $\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$ being the eigenvalues of
$W^{+}$. Let $M_{W^{+}}\subset M$ be the open dense set where the number of
distinct eigenvalues of $W^{+}$ is locally constant. In $M_{W^{+}}$, the
pointwise formula (3.2) is valid locally in the sense that the mutually
orthogonal sections $\omega_{1},\omega_{2},\omega_{3}$ of $\Lambda^{+}$
satisfying (3.1) and the functions $\lambda_{1},\lambda_{2},\lambda_{3}$ may
be assumed differentiable in a neighborhood of any $p\in M_{W^{+}}$.

Since $\Lambda^{+}$ is invariant under parallel transport, in a neighborhood
of $p\in M_{W^{+}}$, there exist one-forms $a,b$, and $c$ defined near $p$,
such that we have (3.2) and

$\displaystyle\nabla\omega_{1}$
$\displaystyle=a\otimes\omega_{2}-c\otimes\omega_{3},$
$\displaystyle\nabla\omega_{2}$
$\displaystyle=b\otimes\omega_{3}-a\otimes\omega_{1},$
$\displaystyle\nabla\omega_{3}$
$\displaystyle=c\otimes\omega_{1}-b\otimes\omega_{2}.$

It was shown by Derdzinski [Der83] that if $\delta W^{+}=0$, then in a
neighborhood of $p\in M_{W^{+}}$, we have

$\displaystyle\nabla\lambda_{1}$ $\displaystyle=$
$\displaystyle(\lambda_{2}-\lambda_{1})(\iota_{a^{\\#}}\omega_{3})^{\\#}+(\lambda_{3}-\lambda_{1})(\iota_{c^{\\#}}\omega_{2})^{\\#},$
$\displaystyle\nabla\lambda_{2}$ $\displaystyle=$
$\displaystyle(\lambda_{1}-\lambda_{2})(\iota_{a^{\\#}}\omega_{3})^{\\#}+(\lambda_{3}-\lambda_{2})(\iota_{b^{\\#}}\omega_{1})^{\\#},$
$\displaystyle\nabla\lambda_{3}$ $\displaystyle=$
$\displaystyle(\lambda_{1}-\lambda_{3})(\iota_{c^{\\#}}\omega_{2})^{\\#}+(\lambda_{2}-\lambda_{3})(\iota_{b^{\\#}}\omega_{1})^{\\#},$

and

$\displaystyle\Delta\lambda_{1}$ $\displaystyle=$ $\displaystyle
2(\lambda_{1}-\lambda_{2})|(\iota_{a^{\\#}}\omega_{3})^{\\#}|^{2}+2(\lambda_{1}-\lambda_{3})|(\iota_{c^{\\#}}\omega_{2})^{\\#}|^{2}$
$\displaystyle+\tfrac{S}{2}\lambda_{1}-2\lambda_{1}^{2}-4\lambda_{2}\lambda_{3},$
$\displaystyle\Delta\lambda_{2}$ $\displaystyle=$ $\displaystyle
2(\lambda_{2}-\lambda_{1})|(\iota_{a^{\\#}}\omega_{3})^{\\#}|^{2}+2(\lambda_{2}-\lambda_{3})|(\iota_{b^{\\#}}\omega_{1})^{\\#}|^{2}$
$\displaystyle+\tfrac{S}{2}\lambda_{2}-2\lambda_{2}^{2}-4\lambda_{1}\lambda_{3},$
$\displaystyle\Delta\lambda_{3}$ $\displaystyle=$ $\displaystyle
2(\lambda_{3}-\lambda_{1})|(\iota_{c^{\\#}}\omega_{2})^{\\#}|^{2}+2(\lambda_{3}-\lambda_{2})|(\iota_{b^{\\#}}\omega_{1})^{\\#}|^{2}$
$\displaystyle+\tfrac{S}{2}\lambda_{3}-2\lambda_{3}^{2}-4\lambda_{1}\lambda_{2},$

where $\iota$ is the interior product and $\\#$ is the sharp operator.

It follows from the assumption $\lambda_{2}\geq-S/12$ that in the set
$M_{W^{+}}$, we have

$\displaystyle\Delta(\lambda_{3}-\lambda_{1})$ $\displaystyle=$ $\displaystyle
6(\lambda_{3}-\lambda_{1})\left(\lambda_{2}+\tfrac{S}{12}\right)+4(\lambda_{3}-\lambda_{1})|(\iota_{c^{\\#}}\omega_{2})^{\\#}|^{2}$
$\displaystyle+2(\lambda_{3}-\lambda_{2})|(\iota_{b^{\\#}}\omega_{1})^{\\#}|^{2}+2(\lambda_{2}-\lambda_{1})|(\iota_{a^{\\#}}\omega_{3})^{\\#}|^{2}$
$\displaystyle\geq$ $\displaystyle
4(\lambda_{3}-\lambda_{1})|(\iota_{c^{\\#}}\omega_{2})^{\\#}|^{2}+2(\lambda_{3}-\lambda_{2})|(\iota_{b^{\\#}}\omega_{1})^{\\#}|^{2}$
$\displaystyle+2(\lambda_{2}-\lambda_{1})|(\iota_{a^{\\#}}\omega_{3})^{\\#}|^{2}.$

Since $\lambda_{3}-\lambda_{1}$ is a nonnegative constant on $M$ by
Proposition 2.1, we conclude that

$(\lambda_{3}-\lambda_{1})|(\iota_{c^{\\#}}\omega_{2})^{\\#}|^{2}=(\lambda_{3}-\lambda_{2})|(\iota_{b^{\\#}}\omega_{1})^{\\#}|^{2}=(\lambda_{2}-\lambda_{1})|(\iota_{a^{\\#}}\omega_{3})^{\\#}|^{2}=0,$

in a neighborhood of $p\in M_{W^{+}}$. This implies that

$\nabla\lambda_{1}=\nabla\lambda_{3}=0$

in that neighborhood of $p\in M_{W^{+}}$. Therefore, $\lambda_{1}$ and
$\lambda_{3}$ are locally constant on $M_{W^{+}}$. Since $M_{W^{+}}$ is open
and dense in $M$, and since $\lambda_{1}$ and $\lambda_{3}$ are locally
Lipschitz functions, we conclude that $\lambda_{1}$ and $\lambda_{3}$ are
global constant functions. Moreover $\lambda_{2}=-\lambda_{1}-\lambda_{3}$ is
also a constant.

If $M$ is not anti-self-dual, then we have $\lambda_{2}=-\frac{S}{12}$ by
Proposition 2.1. $S$ must be a constant since $\lambda_{2}$ is so.

$\square$

###### Proof of Theorem 1.8.

If $M$ is anti-self-dual, then we have $S\geq 0$ in view of
$0=\lambda_{2}\geq-\frac{S}{12}$. If $M$ is not anti-self-dual, then $M$ has
constant scalar curvature by Proposition 3.1. In the latter case, Proposition
2.2 is valid and the proof of Kählerity is the same as the proof of Theorem
1.1 given in Section 2. $\square$

It is clear that Corollary 1.9 and Corollary 1.10 can be observed from the
proof of Theorem 1.8.

## Acknowledgment

The authors would like to thank Professor Jiaping Wang for some helpful
discussions related to this work.

## References

  * [Bes08] Arthur L. Besse.  Einstein manifolds.  Classics in Mathematics. Springer-Verlag, Berlin, 2008.  Reprint of the 1987 edition.
  * [Cal58] E. Calabi.  An extension of E. Hopf’s maximum principle with an application to Riemannian geometry.  Duke Math. J., 25:45–56, 1958.
  * [CLN06] Bennett Chow, Peng Lu, and Lei Ni.  Hamilton’s Ricci flow, volume 77 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI; Science Press Beijing, New York, 2006.
  * [Der83] Andrzej Derdziński.  Self-dual Kähler manifolds and Einstein manifolds of dimension four.  Compositio Math., 49(3):405–433, 1983.
  * [FKP14] Joel Fine, Kirill Krasnov, and Dmitri Panov.  A gauge theoretic approach to Einstein 4-manifolds.  New York J. Math., 20:293–323, 2014.
  * [GL99] Matthew J. Gursky and Claude Lebrun.  On Einstein manifolds of positive sectional curvature.  Ann. Global Anal. Geom., 17(4):315–328, 1999.
  * [Ham86] Richard S. Hamilton.  Four-manifolds with positive curvature operator.  J. Differential Geom., 24(2):153–179, 1986.
  * [Hit74] Nigel Hitchin.  Compact four-dimensional Einstein manifolds.  J. Differential Geometry, 9:435–441, 1974.
  * [LeB15] Claude LeBrun.  Einstein metrics, harmonic forms, and symplectic four-manifolds.  Ann. Global Anal. Geom., 48(1):75–85, 2015.
  * [LeB21] Claude LeBrun.  Einstein manifolds, self-dual Weyl curvature, and conformally Kähler geometry.  Math. Res. Lett., 28(1):127–144, 2021.
  * [MW93] Mario J. Micallef and McKenzie Y. Wang.  Metrics with nonnegative isotropic curvature.  Duke Math. J., 72(3):649–672, 1993.
  * [Pol92] A. Polombo.  De nouvelles formules de Weitzenböck pour des endomorphismes harmoniques. Applications géométriques.  Ann. Sci. École Norm. Sup. (4), 25(4):393–428, 1992.
  * [RS16] Thomas Richard and Harish Seshadri.  Positive isotropic curvature and self-duality in dimension 4.  Manuscripta Math., 149(3-4):443–457, 2016.
  * [Wu17] Peng Wu.  A Weitzenböck formula for canonical metrics on four-manifolds.  Trans. Amer. Math. Soc., 369(2):1079–1096, 2017.
  * [Wu21] Peng Wu.  Einstein four-manifolds with self-dual Weyl curvature of nonnegative determinant.  Int. Math. Res. Not. IMRN, (2):1043–1054, 2021.