3cm3cm3cm3cm

# Subset SSD for enhanced indexation with sector constraints

Cristiano Arbex Valle1 Cristiano Arbex Valle is funded by FAPEMIG grant
APQ-01267-18. John E Beasley2

###### Abstract

In this paper we apply second order stochastic dominance (SSD) to the problem
of enhanced indexation with asset subset (sector) constraints. The problem we
consider is how to construct a portfolio that is designed to outperform a
given market index whilst having regard to the proportion of the portfolio
invested in distinct market sectors.

In our approach, subset SSD, the portfolio associated with each sector is
treated in a SSD manner. In other words in subset SSD we actively try to find
sector portfolios that SSD dominate their respective sector indices. However
the proportion of the overall portfolio invested in each sector is not pre-
specified, rather it is decided via optimisation.

Computational results are given for our approach as applied to the S&P 500
over the period $29^{\text{th}}$ August 2018 to $29^{\text{th}}$ December
2023. This period, over 5 years, includes the Covid pandemic, which had a
significant effect on stock prices. Our results indicate that the scaled
version of our subset SSD approach significantly outperforms the S&P 500 over
the period considered. Our approach also outperforms the standard SSD based
approach to the problem.

1Departamento de Ciência da Computação,
Universidade Federal de Minas Gerais,
Belo Horizonte, MG 31270-010, Brasil
<EMAIL_ADDRESS>
2Brunel University
Mathematical Sciences, UK
<EMAIL_ADDRESS>

Keywords: enhanced indexation, finance, optimisation, portfolio optimisation,
second order stochastic dominance

##  1 Introduction

In this paper we consider the problem of enhanced indexation with asset subset
(sector) constraints. In this problem we aim to outperform a given market
index whilst having regard to the proportion of the portfolio invested in
distinct market sectors. We apply second order stochastic dominance (SSD) to
the problem. Computational results are given for our approach as applied to
the S&P 500.

The structure of this paper is as follows. In Section 2 we review the relevant
literature as to second order stochastic dominance. In Section 3 we present
SSD from a mathematical viewpoint together with discussion of the standard
cutting plane procedure associated with its resolution. In Section 4 we
present our subset SSD approach when we have sector (asset subset) constraints
present that constrain investment in a number of different subsets of assets.
In Section 5 we present computational results obtained when our subset SSD
approach is applied to the S&P 500. In Section 6 we present our conclusions.

We believe that the contribution to the literature of this paper is:

  * •

to present a new approach, _subset SSD_ , for the problem of enhanced
indexation with asset subset (sector) constraints

  * •

to demonstrate computationally, using data that we make publicly available,
that our subset SSD approach significantly outperforms both the S&P 500 and
the standard SSD approach to the problem

##  2 Literature review

The importance of stochastic dominance (SD) within financial portfolio
selection has been recognised for decades (Hadar and Russell, 1969; Bawa,
1975; Levy, 1992). For two random variables $X$ and $Y$ it is well known that
$X$ dominates $Y$ under first-order stochastic dominance (FSD,
$X\succeq_{{}_{FSD}}Y$) if and only if it is preferrable over any monotonic
increasing utility function. Likewise, $X$ dominates $Y$ under second-order
stochastic dominance (SSD, $X\succeq_{{}_{SSD}}Y$) if and only if it is
preferrable over any increasing and strictly concave (risk-averse) utility
function (Whitmore and Findlay, 1978).

For many years however SD was primarily a theoretical framework in terms of
financial portfolio optimisation. This was due to the perceived computational
difficulties associated with finding SD-efficient portfolios. In the past
twenty years, however, there has been a shift towards applying SD (especially
SSD) principles in practice, with several optimisation approaches having been
proposed for finding portfolios that are either SSD-efficient (with regards to
a specified set of feasible portfolios) or SSD-dominating (with regards to a
benchmark).

Ogryczak and Ruszczynski (2002) identified several risk measures that can be
employed in mean-risk ($\mu_{x},r_{X}$) decision models that are consistent
with the SSD relation in the sense that $X\succeq_{{}_{SSD}}Y$ implies that
$\mu_{X}\geq\mu_{Y}$ and $r_{X}\leq r_{Y}$. These measures include tail value-
at-risk, tail Gini mean difference and weighted mean deviation from a
quantile. The authors presented stochastic linear programming formulations for
these models whose optimal solutions are guaranteed to be SSD-efficient.

Kuosmanen (2004, 2001) developed the first SSD efficiency tests based on
mathematical programming. Their formulation finds, if it exists, the portfolio
with the highest in-sample mean that dominates a benchmark in the SSD sense.
Post (2003) developed linear programming models for testing if a given
portfolio is SSD-efficient with respect to all possible portfolios given a set
of assets.

Dentcheva and Ruszczynski (2006, 2003) first combine the available assets to
produce a reference (or benchmark) distribution, and then compute a portfolio
which SSD-dominates the benchmark. They used the lower partial moment of order
one to develop the SSD ranking concerning the benchmark portfolio. Their work
has been the basis of several later papers in literature, as referenced below.

Roman et al. (2006) introduced a multi-objective optimisation model to find a
portfolio that achieves SSD dominance over a benchmark. If no such portfolio
exists they find the portfolio whose return distribution comes closest to the
benchmark. They showed that SSD efficiency does not necessarily make a return
distribution desirable, as demonstrated by the optimal portfolio with regards
to maximum expected return (which is SSD-efficient). They emphasised the
crucial role played by a carefully selected benchmark in the process.

Luedtke (2008) presented a model that generalises that of Kuosmanen (2004)
which includes FSD constraints based on a cutting-plane formulation for
problems with integrated chance constraints. Their model involves integer
variables, but relaxing integrality yields a formulation with SSD constraints.
Their objective is to maximise expected portfolio return.

Fábián et al. (2011a, b) introduced a cutting plane reformulation of Roman et
al. (2006) which generalises Dentcheva and Ruszczynski (2006). The authors
replaced the multi-objective nature of the problem by maximising the minimum
value in the SSD relation with regards to a benchmark. Roman et al. (2013)
applied the SSD cutting plane formulation in an enhanced indexation setting.
Valle et al. (2017) added exogenous constraints and reformulated the problem
as an integer linear program, for which a branch-and-cut algorithm was
developed.

Kopa and Post (2015); Post and Kopa (2013) introduced a more generalised
efficiency test which allows for unequal probabilities and higher orders. In
the case of inefficiency their dual model finds a dominating portfolio. If the
portfolio being tested is a benchmark, this dual model can be seen as
equivalent to a model for enhanced indexation.

The set of SSD efficient portfolios is generally very large, and investors
need to decide how to select a portfolio in which to invest from within this
set. The formulation from Post and Kopa (2013) may be used to find different
SSD-efficient portfolios depending on how some parameters are specified.
Hodder et al. (2015) proposed ways to assign values to these parameters with
the goal of helping investors select a single portfolio out of the efficient
set.

Bruni et al. (2017, 2012) developed an alternative approach for SD-based
enhanced indexation. They proposed a criterion called “cumulative zero-order
stochastic $\epsilon$-dominance” (CZS$\epsilon$D). Zero-order SD happens when
all returns from a given portfolio are superior to all returns from an
alternative portfolio. The authors attempt to minimise underperformance by
adding an exponential number of constraints related to the CZS$\epsilon$D
criterion, where $\epsilon$ is the maximum underperformance allowed. The
separation algorithm they use is equivalent to optimising conditional value-
at-risk via linear programming.

Sharma et al. (2017) introduced a relaxed-SSD formulation for enhanced
indexation. The SSD constraints are relaxed by adding under/overachievement
where SSD violation is controlled by setting an appropriate upper bound
related to the total underachievement. The concept of relaxed-SSD was first
introduced by Lizyayev and Ruszczyński (2012).

Sharma and Mehra (2017) proposed a SSD-based approach for producing sector
portfolios. For each sector, their model seeks a SSD portfolio that dominates
the corresponding sector index, whilst focusing on a number of financial
ratios when making sector portfolio decisions. These sector portfolios are
then combined using another model that optimises their mean return subject to
being (if possible) SSD-dominating with respect to the main market index. If
SSD dominance cannot be achieved, either in relation to a sector, or in
relation to the main market index, they relax the dominance constraints in
their models.

Liu et al. (2021) showed that FSD and SSD may not be sufficient to
discriminate between multiple dominating portfolios with regards to a
benchmark. They proposed a new criterion called Interval-based SD (ISD) in
which different SD orders are applied to different parts of the support of the
return distribution. They present a reformulation of Dentcheva and Ruszczynski
(2006) that maximises portfolio return subject to ISD constraints.

Sehgal and Mehra (2021) presented a robust version of the SSD-formulation of
Dentcheva and Ruszczynski (2006). Robustness is introduced by varying asset
returns, and the model is developed as the deterministic equivalent of a
stochastic programming formulation. Goel and Sharma (2021) also generalised
Dentcheva and Ruszczynski (2006) by considering the “utility improvement” in
portfolio returns instead of the returns themselves. The authors proposed
replacing the portfolio and benchmark returns by their respective deviations
in the SSD constraints.

Malavasi et al. (2021) compared the performance of SSD portfolios with
efficient portfolios derived using the standard mean-variance approach of
Markowitz (1952). They also focused on the performance of the global minimum
variance portfolio as compared with portfolios that are stochastically
dominant to this minimum variance portfolio.

Cesarone et al. (2023) compared the formulations of Roman et al. (2013) and
Kopa and Post (2015) with skewed benchmarks obtained by using the reshaping
method of Valle et al. (2017). They found that SSD portfolios that dominate
the skewed benchmark generally perform better out-of-sample.

Liesio et al. (2023) considered the problem of generating an efficient
frontier using stochastic dominance. They presented an approach based on
Pareto optimal solutions of a multiple objective optimisation problem.

Cesarone and Puerto (2024) presented an alternative to Roman et al. (2013)
where, instead of maximising the minimum value of the SSD relation, the
authors proposed a model that optimises the ordered weighted average of a
predefined number of tails.

##  3 Cutting plane based SSD formulation

Based on a reformulation of the conditional value-at-risk minimisation problem
given by Künzi-Bay and Mayer (2006), Fábián et al. (2011a) proposed a novel
cutting plane formulation of the SSD problem, one whose objective is to
maximise the minimum value in the SSD relationship between the portfolio and a
given benchmark (e.g. a market index or some reference distribution). Roman et
al. (2013) then employed the formulation for enhanced indexation. In this
section we outline their approach. Let

  * •

$N$ be number of assets available for investment

  * •

$S$ be number of scenarios, where the scenarios are assumed to be equiprobable

  * •

$r_{is}$ be the return of asset $i$ in scenario $s$

  * •

$r^{I}_{s}$ be the benchmark return in scenario $s$

  * •

$R_{s}^{P}$ be the return associated with a given asset portfolio $P$ in
scenario $s$

  * •

$\text{Tail}^{L}_{\frac{\alpha}{S}}(P)$ be the unconditional expectation of
the smallest $\alpha$ outcomes in $[R_{s}^{P}~{}|~{}s=1,\ldots,S]$, so the
left tail of the portfolio return distribution

For a portfolio $P$ with asset weights $[w_{i}]$ and hence return
$R_{s}^{P}=\sum_{i=1}^{N}r_{is}w_{i}$ in scenario $s$ we, as Fábián et al.
(2011a) albeit with slightly different notation, define
$\text{Tail}^{L}_{\frac{s}{S}}(P)$ using

$\text{Tail}^{L}_{\frac{s}{S}}(P)=\frac{1}{S}\text{(sum of the $s$ smallest
portfolio returns in $[R_{1}^{P},R_{2}^{P},\ldots,R_{S}^{P}]$)}$ (1)

Here $\text{Tail}^{L}_{\frac{s}{S}}(P)$ is the left tail of the cumulative
return distribution associated with $[R_{1}^{P},R_{2}^{P},\ldots,R_{S}^{P}]$
weighted by the constant $(1/S)$ factor.

Let $I$ be some index portfolio which we would (ideally) like to outperform.
The index portfolio has known return $R_{s}^{I}$ in scenario
$s,~{}s=1,\ldots,S$. Let
$\hat{\tau}_{s}=\text{Tail}^{L}_{\frac{s}{S}}(I)~{}s=1,\ldots,S$. Clearly we
would like the tails of the chosen portfolio to improve on the index portfolio
tails, so define the tail differences $\mathcal{V}_{s}$ between the chosen
portfolio and the index portfolio using

$\mathcal{V}_{s}=\text{Tail}^{L}_{\frac{s}{S}}(P)-\hat{\tau}_{s}~{}~{}~{}s=1,\ldots,S$
(2)

If $\mathcal{V}_{s}\geq 0~{}s=1,\ldots,S$ then the portfolio is second order
stochastic dominant to the index portfolio.

Now it is trivial to observe that the sum of the $s$ smallest portfolio
returns in the $S$ scenarios can be found by considering all subsets
$\mathcal{J}$ of the $S$ scenarios of cardinality $s$. In other words

$\text{Tail}^{L}_{\frac{s}{S}}(P)=\frac{1}{S}\min\left[\sum_{j\in\mathcal{J}}\sum_{i=1}^{N}r_{ij}w_{i}~{}|~{}\mathcal{J}\subseteq\\{1,...,S\\},|\mathcal{J}|=s\right]$
(3)

If we are choosing $s$ scenarios from the $S$ scenarios then there are
$\frac{S!}{s!(S-s)!}$ subsets $\mathcal{J}$ that need to be considered. So
Equation (3) defines the $s$ smallest portfolio returns in the $S$ scenarios
using a combinatorial number of constraints.

Now to make use of the combinatorial definition of
$\text{Tail}^{L}_{\frac{s}{S}}(P)$ let $\mathcal{V}$ be the minimum value of
$[\mathcal{V}_{s}~{}|~{}s=1,\ldots,S]$. Then a suitable optimisation program
to decide the portfolio of assets that should be held is to

$\mbox{maximise}~{}\mathcal{V}$ (4)

subject to

$\mathcal{V}_{s}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i=1}^{N}r_{ij}w_{i}-\hat{\tau}_{s}~{}~{}~{}~{}\forall\mathcal{J}\subseteq\\{1,...,S\\},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S$
(5) $\mathcal{V}\leq\mathcal{V}_{s}~{}~{}~{}~{}s=1,\ldots,S$ (6)
$\sum_{i=1}^{N}w_{i}=1$ (7) $0\leq w_{i}\leq 1~{}~{}~{}~{}i=1,\ldots,N$ (8)
$\mathcal{V}\in\mathbb{R}$ (9)
$\mathcal{V}_{s}\in\mathbb{R}~{}~{}~{}s=1,\ldots,S$ (10)

Equation (4), in conjunction with Equation (6), maximises the minimum tail
difference. Equation (5) is the standard SSD combinatorial definition of the
tail differences. Equation (7) ensures that all of our wealth is invested in
assets. Equation (8) is the non-negativity constraint (so no short-selling).
Equation (9) ensures that $\mathcal{V}$ can be positive or negative whilst
Equation (10) ensures that the tail differences $\mathcal{V}_{s}$ can be
positive or negative.

Equations (4)-(10) above is a portfolio choice optimisation program with
explicit consideration of tails. If the objective function has a non-negative
optimal value then the associated portfolio is second order stochastic
dominant with respect to the index.

###  3.1 Cutting plane resolution

We can adopt a cutting plane resolution procedure for the portfolio
optimisation program Equations (4)-(10) above. This has been given previously
(albeit in a slightly different form) by Fábián et al. (2011a).

First define an initial scenario set $\mathcal{J^{*}}$ where there is at least
one set of cardinality $s$, for all values of $s=1,\ldots,S$, in
$\mathcal{J^{*}}$ and amend Equation (5) to

$\mathcal{V}_{s}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i=1}^{N}r_{ij}w_{i}-\hat{\tau}_{s}~{}~{}~{}~{}\forall\mathcal{J}\in\mathcal{J^{*}},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S$
(11)

  1. 1.

Solve the amended optimisation program, optimise Equation (4) subject to
Equations (6)-(11).

  2. 2.

Consider each value of $s$ ($s=1,\ldots,S$) in turn and if in the solution to
the amended optimisation program

$\mathcal{V}_{s}>\frac{1}{S}\text{(sum of the $s$ smallest portfolio returns
over the $S$ scenarios)}-\hat{\tau}_{s}$ (12)

then add the scenario set associated with these $s$ smallest portfolio returns
to $\mathcal{J^{*}}$. Here the scenario set that is added constitutes a valid
cut associated with Equation (5) that is violated by the current solution.

  3. 3.

If scenarios sets have been added to $\mathcal{J^{*}}$ go to Step (1), else
terminate.

Upon termination at Step (3) above we will have a set of values satisfying all
of the constraints in the amended optimisation program. It remains to prove
that we have solved the original (unamended) optimisation program to
optimality. Here the only difference between the original optimisation program
and the amended optimisation program is the replacement of Equation (5) by
Equation (11).

Consider a particular value of $s$. Since we have terminated no cuts of the
form shown in Equation (12) can be added, in other words we must have

$\mathcal{V}_{s}\leq\frac{1}{S}\text{(sum of the $s$ smallest portfolio
returns over the $S$ scenarios)}-\hat{\tau}_{s}$ (13)

But the term (sum of the $s$ smallest portfolio returns over the $S$
scenarios) corresponds to
$\min[\sum_{j\in\mathcal{J}}\sum_{i=1}^{N}r_{ij}w_{i}~{}|~{}\mathcal{J}\subseteq\\{1,...,S\\},|\mathcal{J}|=s]$,
since it is the sum of the $s$ smallest portfolio returns. So Equation (13) is
equivalent to

$\mathcal{V}_{s}\leq\frac{1}{S}\min\left[\sum_{j\in\mathcal{J}}\sum_{i=1}^{N}r_{ij}w_{i}~{}|~{}\mathcal{J}\subseteq\\{1,...,S\\},|\mathcal{J}|=s\right]-\hat{\tau}_{s}$
(14)

Equation (14) in turn implies that $\mathcal{V}_{s}$ satisfies Equation (5) in
the original optimisation program. This is because the summation term on the
right-hand side of that equation is over all subsets of cardinality $s$, so
equivalent to the minimisation term in Equation (14). Hence we have found the
optimal solution to the original (unamended) optimisation program.

###  3.2 Scaled tails

One issue with using Equation (4) as an objective is that there may be
multiple distinct portfolios, each of which has the same maximum $\mathcal{V}$
value. However the SSD formulation can be tailored to focus on certain aspects
of the return distribution associated with the portfolio chosen.

With Equation (6) and an objective of maximising $\mathcal{V}$ more importance
is given to $\text{Tail}^{L}_{\frac{s}{S}}(P)$ when $s$ is small. Namely,
$\text{Tail}^{L}_{\frac{s}{S}}(P)$ for $s$ approaching $S$ is given the same
relative importance by Equation (6) as for $s$ close to 1. But since the left
tails are cumulative, for large values of $s$ the most positive portfolio
returns are “diluted” among smaller returns. An unintended consequence of this
is that solving the maximise $\mathcal{V}$ formulation tends to yield
portfolios that have a smaller left tail when compared to benchmark returns
$[\hat{\tau}_{s}~{}|~{}s=1,\ldots,S$], but also a smaller right tail.

As an alternative Fábián et al. (2011b) proposed scaling the tails by
replacing Equation (6) with

$\frac{s}{S}\mathcal{V}\leq\mathcal{V}_{s}~{}~{}~{}~{}s=1,\ldots,S$ (15)

Here the effect of scaling is that more importance is given to the returns in
the right tails of the distribution.

##  4 Subset SSD

Above we have a single set of assets and we seek a portfolio chosen from these
assets that, in a SSD sense, outperforms (if possible) a given asset index. In
this section we generalise this approach to the case where it is possible to
subdivide the entire set of assets into individual subsets, each with
differing characteristics.

We might be interested in different asset subsets for a number of reasons,
e.g. in a given set of index assets it could be that we believe that large
capitalisation assets and low capitalisation assets exhibit different
behaviour. So in our chosen portfolio we might wish to tailor our exposure to
these two different asset subsets differently. Other asset subsets can be
easily envisaged e.g. based on different market sectors, different momentum
characteristics or any other economic metric. In our approach we do not assume
that the asset subsets are disjoint, in other words a single asset can be in
two or more subsets.

_We should be clear here that under the standard SSD approach exposure to
different asset subsets can be included by adding additional constraints to
the SSD formulation, Equations ( 4)-(10), as seen above. However in our
approach EACH individual asset subset portfolio is treated in a SSD manner._
For clarity this standard SSD approach is given in Section 4.2 below.

Suppose that we have $K$ asset subsets where $N^{k}$ are the assets in asset
subset $k$ and $\cup^{K}_{k=1}N^{k}=[1,...,N]$. We need for each asset subset
an underlying index in order to create an appropriate SSD formulation. Such an
index may be publicly available. If not, one can easily be produced using
weights associated with any index that includes these assets.

As an illustration of this suppose that the weight associated with asset $i$
in an appropriate benchmark index is $\Gamma_{i}$, where the index is price
based, so the price $P_{it}$ of asset $i$ at time $t$ contributes to the
index. Then the sub-index for the set $N^{k}$ at time $t$ is given by
$\sum_{i\in N^{k}}\Gamma_{i}P_{it}$, so the index return associated with asset
subset $k$ at time $t$ is $\left[\sum_{i\in N^{k}}\Gamma_{i}P_{it}/\sum_{i\in
N^{k}}\Gamma_{i}P_{it-1}\right]$. Let $I^{k}$ represent the returns on the
index associated with asset subset $k$. Then
$\hat{\tau}^{k}=(\hat{\tau}^{k}_{1},\ldots,\hat{\tau}^{k}_{S})=\big{(}\text{Tail}^{L}_{\frac{1}{S}}I^{k},\ldots,\text{Tail}^{L}_{\frac{S}{S}}I^{k}\big{)}$.

In the SSD formulation below we add a $k$ superscript associated with asset
subset $N^{k}$ to the previous formulation (Equations (4)-(10)). Let
$\mathcal{V}_{s}^{k}$ be the tail difference between the chosen portfolio and
the index portfolio associated with asset subset $k$.

Let $W^{k}\geq 0$ be the proportion of the portfolio invested in subset $k$,
_where this proportion will be decided by the optimisation_. However note here
that, as will be seen below, the decision maker has the flexibility to impose
bounds on $W^{k}$, or indeed to specify the exact value that $W^{k}$ should
take.

Then, drawing on the program given above, Equations (4)-(10), the constraints
of the subset SSD optimisation program are

$\mathcal{V}_{s}^{k}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i\in
N^{k}}r_{ij}w_{i}/W^{k}-\hat{\tau}_{s}^{k}~{}~{}~{}~{}\forall\mathcal{J}\subseteq\\{1,...,S\\},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S,~{}k=1,\ldots,K$
(16) $W^{k}=\sum_{i\in N^{k}}w_{i}~{}~{}~{}~{}k=1,\ldots,K$ (17)
$\delta^{L}_{k}\leq W_{k}\leq\delta^{U}_{k}~{}~{}~{}~{}k=1,\ldots,K$ (18)
$\sum_{i=1}^{N}w_{i}=1$ (19) $0\leq w_{i}\leq 1~{}~{}~{}~{}i=1,\ldots,N$ (20)
$\mathcal{V}_{s}^{k}\in\mathbb{R}~{}~{}~{}~{}s=1,\ldots,S,~{}k=1,\ldots,K.$
(21)

Equation (16) is the tail difference for each subset $k$. In this equation the
summation in the numerator of the first term on the right-hand side of the
inequality is the return from the investment in assets associated with subset
$k$. But unlike Equation (5) above we do not necessarily have that the sum of
the weights (over assets $i\in N^{k}$) will equal one, so we have to scale
this summation by the $W_{k}$ factor before subtracting the
$\hat{\tau}_{s}^{k}$ associated with subset $k$.

Equation (17) defines the subset proportion based on the sum of the
proportions of the total wealth invested in the assets in the subset. Equation
(18) ensures that the proportion of the total investment in subset $k$ lies
between $\delta^{L}_{k}$ and $\delta^{U}_{k}$ where these are the user defined
lower and upper limits on the proportion of the portfolio invested in subset
$k$. Equation (19) ensures that all of our wealth is invested in assets.
Equation (20) is the non-negativity constraint (so no short-selling) and
Equation (21) ensures that the tail differences $\mathcal{V}^{k}_{s}$ can be
positive or negative.

Now assuming that $W^{k}>0~{}k=1,\ldots,K$ (which we can ensure if we wish by
adding constraints $W^{k}\geq\epsilon~{}k=1,\ldots,K$, where $\epsilon{>}0$
and small) we can linearise Equation (16) to

$W^{k}\mathcal{V}_{s}^{k}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i\in
N^{k}}r_{ij}w_{i}–W^{k}\hat{\tau}_{s}^{k}~{}~{}~{}~{}\forall\mathcal{J}\subseteq\\{1,...,S\\},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S,~{}k=1,\ldots,K$
(22)

Here the $W^{k}\mathcal{V}_{s}^{k}$ term is nonlinear, but can be interpreted
as the _proportion weighted tail difference_ associated with set $k$.

Now based on Equation (4) we might be tempted to have an objective function of
the form $\mbox{maximise}~{}\mathcal{V}$ where
$\mathcal{V}\leq\mathcal{V}_{s}^{k}~{}s{=}1,\ldots,S,~{}k{=}1,\ldots,K$ and
$\mathcal{V}\in\mathbb{R}$. Here each tail difference $\mathcal{V}_{s}^{k}$
influences the objective, bounding it from above. However we have _no prior
knowledge of the investment proportion associated with subset $k$_. So for
example if we adopt an objective of this form we might have two subsets with
the same tail difference (as calculated using Equation (16)), so with the same
influence on $\mathcal{V}$, but very different investment proportions.

This seems perverse - surely an investment with a higher proportion should
have more influence with respect to the objective? In other words (somehow)
the investment proportion $W^{k}$ for subset $k$ should ideally be
incorporated, so that the higher the value of $W^{k}$ the more impact subset
$k$ has on the maximisation objective.

It is clear that one way forward is to replace the nonlinear proportion
weighted tail difference term $W^{k}\mathcal{V}_{s}^{k}$ in Equation (22) by a
single term, say $\mathcal{Z}_{s}^{k}\in\mathbb{R}$, and adopt an objective
function of the form

$\mbox{maximise}~{}\mathcal{V}$ (23)

subject to

$\mathcal{Z}_{s}^{k}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i\in
N^{k}}r_{ij}w_{i}–W^{k}\hat{\tau}_{s}^{k}~{}~{}~{}~{}\forall\mathcal{J}\subseteq\\{1,...,S\\},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S,~{}k=1,\ldots,K$
(24)
$\beta\mathcal{V}\leq\mathcal{Z}_{s}^{k}~{}~{}~{}~{}\forall\mathcal{J}\subseteq\\{1,...,S\\},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S,~{}k=1,\ldots,K$
(25) $\mathcal{V}\in\mathbb{R}$ (26)
$\mathcal{Z}_{s}^{k}\in\mathbb{R}~{}~{}~{}~{}s=1,\ldots,S,~{}k=1,\ldots,K.$
(27)

with the other constraints (Equations (17)-(20)) remaining as before. In
Equation (25) $\beta$ is the scaling factor where $\beta{=}1$ for no scaling
and $\beta{=}s/S$ for scaled tails as in Equation (15).

###  4.1 Cutting plane resolution - $\mathcal{Z}_{s}^{k}$

We can adopt what is effectively the same cutting plane resolution procedure
for the portfolio optimisation program as given previously by Fábián et al.
(2011a) and seen above. For completeness here we set out this procedure in
full.

First define an initial scenario set $\mathcal{J^{*}}$ where there is at least
one set of cardinality $s$, for all values of $s=1,\ldots,S$, in
$\mathcal{J^{*}}$ and amend Equation (24) to

$\mathcal{Z}_{s}^{k}\leq\frac{1}{S}\sum_{j\in\mathcal{J}}\sum_{i\in
N^{k}}r_{ij}w_{i}–W^{k}\hat{\tau}_{s}^{k}~{}~{}~{}~{}\forall\mathcal{J}\in\mathcal{J^{*}},~{}|\mathcal{J}|=s,~{}s=1,\ldots,S,~{}k=1,\ldots,K$
(28)

  1. 1.

Solve the amended optimisation program, optimise Equation (23) subject to
Equations (17)-(20),(25)-(28)

  2. 2.

Consider all values of $s$ and $k$ ($s=1,\ldots,S,~{}k=1,\ldots,K$) in turn
and if in the solution to the amended optimisation program

$\mathcal{Z}_{s}^{k}>\frac{1}{S}(\mbox{sum of the $s$ smallest portfolio
returns in subset $k$ over the $S$ scenarios)}–W^{k}\hat{\tau}_{s}^{k}$ (29)

then add the scenario set associated with these $s$ smallest returns to
$\mathcal{J^{*}}$. Here the scenario set that is added constitutes a valid cut
associated with Equation (24) that is violated by the current solution.

  3. 3.

If scenarios sets have been added to $\mathcal{J^{*}}$ go to Step (1), else
terminate.

###  4.2 Standard SSD based approach

Above we have presented our approach where each individual asset subset is
treated in a SSD manner. The standard approach to the problem of how to
construct a portfolio that is designed to outperform a given market index
whilst having regard to the proportion of the portfolio invested in distinct
asset subsets (market sectors) is to add constraints related to asset subsets
to the standard SSD formulation.

In term of the notation given above this approach would correspond to optimise
Equation (4) subject to Equations (5)-(10),(17),(18). Here we have added the
subset constraints, Equations (17),(18), to the standard SSD formulation.

Computational results, presented below, indicate that for the S&P 500 over the
period which we considered, this standard approach is outperformed by our
approach.

##  5 Computational results

We used a dataset associated with the S&P 500, with daily stock prices from
$29^{\text{th}}$ August 2018 until $29^{\text{th}}$ December 2023. This time
period, over 5 years, includes the Covid pandemic, which had a significant
effect on stock prices. Our data has been manually adjusted to account for
survivorship bias - on a given date only assets that were part of the S&P 500
index at that time are available to be selected for investment.

In order to define the scenarios required by SSD we used a lookback approach
that included the most recent 85 daily prices, which then yield 84 in-sample
returns (roughly a quadrimester in business days).

The SSD subsets were defined by the economic sectors to which each asset
belongs. There are 11 different stock market sectors according to the most
commonly used classification system, known as the Global Industry
Classification Standard (GICS). These sectors are communication services,
consumer discretionary, consumer staples, energy, financials, healthcare,
industrials, materials, real estate, technology and utilities. For each
sector, its benchmark consisted of the corresponding time series for the S&P
sector indices111https://www.spglobal.com/spdji/en/index-family/equity/us-
equity/sp-sectors/. Table 1 shows the S&P 500 sector breakdown as of
$9^{\text{th}}$ October 2023 together with the approximate weight of the
sector with regard to the index.

Sector | Approximate weight (%)
---|---
Technology | 26.0
Healthcare | 14.5
Financials | 12.9
Consumer discretionary | 9.9
Industrials | 8.6
Communication services | 8.2
Consumer staples | 7.4
Energy | 4.5
Utilities | 2.9
Materials | 2.6
Real estate | 2.5
Table 1: S&P 500 sector breakdown

All of the data used in this paper is publicly available for the use by other
researchers at:

https://github.com/cristianoarbex/subsetSSDData/

We used CPLEX Optimizer 22.1.0 (2023) as the linear and integer programming
solver, with default options. Our backtesting tool is developed in Python and
all optimisation models are developed in C++. We ran all experiments on an
Intel(R) Core(TM) i7-3770 CPU @ 3.90GHz with 8 cores, 8GB RAM and with Ubuntu
22.04.3 LTS as the operating system.

###  5.1 Out-of-sample performance

In this section we evaluate the performance of our subset SSD approach when
compared to both the S&P 500 and the standard SSD approach with sector
constraints, which was outlined above in Section 4.2.

As mentioned above we used an in-sample period of 85 days. We conducted
periodic rebalancing every 21 days (roughly one month in business days). To
illustrate our approach our first in-sample period of 85 days runs from
$29^{\text{th}}$ August 2018 until $31^{\text{st}}$ December 2018. So using
this in-sample period (with 84 return values for each asset) we choose a
portfolio (using a SSD strategy) on $31^{\text{st}}$ December 2018, evaluate
its performance out-of-sample for the next 20 business days so from
$31^{\text{st}}$ December 2018 to $31^{\text{st}}$ January 2019, then repeat
the process until the data is exhausted. In total this involved 60 out-of-
sample periods for which we then have a single out-of-sample time series of
returns. For simplicity we assume no transaction costs.

We evaluated four different strategies. The unscaled and scaled versions of
our subset SSD approach and, equivalently, the unscaled and scaled versions of
the standard SSD approach with sector constraints.

In order to define sector bounds, for a given sector $k$ we take its exposure
from Table 1 as $\delta_{k}$ and define an interval $\Delta=0.05$ such
$\delta_{k}^{L}=(1-\Delta)\delta_{k}$ and
$\delta_{k}^{U}=(1+\Delta)\delta_{k}$, where $\delta_{k}^{L}$ and
$\delta_{k}^{U}$ limit exposure to any particular sector, as in Equation (18).
These bounds apply to both subset SSD and standard SSD. This choice ensures
that the portfolios chosen under both subset and standard SSD have similar
exposure to S&P 500 sectors, whilst, at the same time, giving some leeway to
the SSD optimiser in its choice of portfolio.

Figures 2 and 2 show graphically the cumulative returns during the out-of-
sample period for all four strategies and the S&P 500. For easier
visualisation, we show these results separately, with the subset SSD results
in Figure 2 and the standard SSD results in Figure 2. Both figures use exactly
the same scale.

Figure 1: Cumulative out-of-sample returns for both the unscaled and scaled
versions of the subset SSD formulation with sector constraints

Figure 2: Cumulative out-of-sample returns for both the unscaled and scaled
versions of the standard SSD formulation with sector constraints

Considering these figures the effect of the Covid pandemic can be clearly
seen, with a dramatic fall in cumulative returns for the S&P 500 in the first
half of 2020. It is clear from these figures that the scaled version of our
subset SSD approach significantly outperforms the S&P 500 over the time period
considered.

In order to gain some numeric insight into the performance of the four
strategies as seen in Figure 2 and Figure 2 we show some selected comparative
statistics in Table 2. These are calculated from the out-of-sample returns for
the four strategies, and correspondingly for the S&P 500 index.

Let $Q$ be a series of $0,\ldots,T$ daily portfolio values, where $Q_{t}$ is
the value of the given portfolio on day $t$. In Table 2 FV stands for the
final portfolio value, assuming a starting amount of $1, and is calculated as
$Q_{T}/Q_{0}$. CAGR stands for Capital Annualised Growth Rate and as a
percentage is calculated as
$100\left(\left(\frac{Q_{T}}{Q_{0}}\right)^{\frac{1}{Y}}-1\right)$, where
$Y=T/252$ is an approximation for the number of years in the out-of-sample
period. Column Vol represents the annualised sample standard deviation of the
out-of-sample returns. Sharpe and Sortino are the annualised Sharpe and
Sortino ratios respectively, where for their calculation we use the CBOE
10-year treasury notes (symbol TNX) as the risk-free rate. MDD represents the
maximum drawdown and as a percentage is calculated as
$\max\left(0,100\max_{0\leq t<u\leq T}\frac{Q_{t}-Q_{u}}{Q_{t}}\right)$.

Strategies | FV | CAGR | Vol | Sharpe | Sortino | MDD
---|---|---|---|---|---|---
Subset SSD (scaled) | 2.28 | 17.92 | 18.57 | 0.86 | 1.22 | 29.05
Subset SSD (unscaled) | 1.57 | 9.40 | 19.53 | 0.44 | 0.60 | 36.30
Standard SSD (scaled) | 1.80 | 12.48 | 20.89 | 0.56 | 0.77 | 33.93
Standard SSD (unscaled) | 1.58 | 9.55 | 17.15 | 0.49 | 0.67 | 28.08
S&P 500 | 1.90 | 13.74 | 21.31 | 0.61 | 0.84 | 33.92
Table 2: Comparative out-of-sample statistics

With regard to the scaling of tails, Fábián et al. (2011b); Roman et al.
(2013); Valle et al. (2017) all concluded that scaled SSD tends to achieve
superior out-of-sample returns, but not necessarily superior risk, when
compared to unscaled SSD. The reason for this is that by scaling the tails
more importance is given to the returns in the right tails of the
distribution. Here we observe the same behaviour, with the scaled versions of
both standard and subset SSD outperforming their unscaled versions in terms of
performance (FV, CAGR). The gain in absolute performance also translates to
better risk-adjusted performance (Sharpe, Sortino). As can be seen from Table
2 the unscaled formulations both show inferior performance when compared to
the S&P 500.

With regards to the scaled formulations, subset SSD performed considerably
better than standard SSD with sector constraints. We would remind the reader
here that the main difference between the two approaches is that with subset
SSD we actively try to find sector portfolios that SSD dominate their
respective sector indices, as opposed to standard SSD where there is no
attempt to ensure this.

Subset SSD achieved better returns (in terms of FV and CAGR) and better risk
(in terms of Vol and MDD) and therefore much improved risk-adjusted
performance (Sharpe, Sortino) as compared with standard SSD and too as
compared with the S&P 500. Despite the Covid drop in 2020, during the entire
period considered the S&P 500 had a strong positive performance (almost
doubling in value). However subset SSD was able not only to outperform the S&P
500 in terms of return, but also in terms of risk.

Despite the potentially exponential number of constraints involved in the
cutting plane procedures for SSD solution our experience has been that the
computational effort required to solve each portfolio rebalance to optimality
was negligible. In our experiments a total of 60 rebalances were needed. For
the scaled subset SSD formulation the average computational time per rebalance
was 0.58s, with a maximum of 1.86s and a minimum of 0.13s (median 0.54s),
while for the other strategies the average computational time was between 0.3s
and 0.35s and no rebalance required more than a second.

Figure 3 shows the exposure per sector for scaled subset SSD. The figure shows
comparatively little variation per sector, as expected, since the strategies
are limited by sector bounds to be within $\Delta$, here 5%, of the sector
weightings in the S&P 500.

Figure 3: Out-of-sample exposure per sector, scaled subset SSD

###  5.2 Varying sector bounds

To investigate the performance of our subset SSD approach when we varied
sector bounds we performed ten different experiments. As above, in order to
define sector bounds for a given sector $k$ we take its exposure from Table 1
as $\delta_{k}$. Using $\Delta$ we have $\delta_{k}^{L}=(1-\Delta)\delta_{k}$
and $\delta_{k}^{U}=(1+\Delta)\delta_{k}$, where (as before) $\delta_{k}^{L}$
and $\delta_{k}^{U}$ limit exposure to any particular sector, as in Equation
(18).

We evaluated the out-of-sample performance of both scaled subset SSD and
scaled standard SSD for $\Delta=(0.01,0.02,\dots,0.10)$. The results can be
seen in Table 3. In this table we have, for example for FV and scaled subset
SSD, that over the ten values of $\Delta$ considered, the mean FV value was
2.18, the median FV value was 2.13, the minimum FV value was 1.97 and the
maximum FV value was 2.37.

It is clear from Table 3 that, for the data we considered, scaled subset SSD
is superior to scaled standard SSD. For the four performance measures where
high values are better (so FV, CAGR, Sharpe and Sortino) the _minimum_ values
for these measures for scaled subset SSD exceed the _maximum_ values for these
measures for scaled standard SSD. For the two performance measures where low
values are better (so Vol and MDD) the _maximum_ values for these measures for
scaled subset SSD are below the _minimum_ values for these measures for scaled
standard SSD. _In other words with regard to all six performance measures
scaled subset SSD dominates scaled standard SSD._

In a similar fashion for the four performance measures where high values are
better (so FV, CAGR, Sharpe and Sortino) the minimum values for these measures
for scaled subset SSD exceed the values associated with the S&P 500. For the
two performance measures where low values are better (so Vol and MDD) the
maximum values for these measures for scaled subset SSD are below the values
associated with the S&P 500. _In other words with regard to all six
performance measures scaled subset SSD dominates the S &P 500._

Stats | Subset SSD (scaled) | Standard SSD (scaled) | S&P 500
---|---|---|---
Mean | Median | Min | Max | Mean | Median | Min | Max |
FV | 2.18 | 2.13 | 1.97 | 2.37 | 1.80 | 1.80 | 1.80 | 1.81 | 1.90
CAGR | 16.88 | 16.32 | 14.52 | 18.90 | 12.51 | 12.51 | 12.44 | 12.59 | 13.74
Vol | 18.55 | 18.57 | 18.31 | 18.74 | 20.91 | 20.91 | 20.82 | 21.02 | 21.31
Sharpe | 0.81 | 0.79 | 0.70 | 0.91 | 0.56 | 0.56 | 0.56 | 0.57 | 0.61
Sortino | 1.15 | 1.12 | 0.97 | 1.29 | 0.77 | 0.77 | 0.77 | 0.78 | 0.84
MDD | 29.29 | 29.05 | 28.35 | 30.87 | 33.93 | 34.02 | 33.54 | 34.27 | 33.92
Table 3: Summary statistics for the scaled formulations when
$\Delta=(0.01,0.02,\ldots,0.10)$

##  6 Conclusions

In this paper we have considered the problem of how to construct a portfolio
that is designed to outperform a given market index, whilst having regard to
the proportion of the portfolio invested in distinct market sectors.

We presented a new approach, subset SSD, for the problem. In our approach
portfolios associated with each sector are treated in a SSD manner so that we
actively try to find sector portfolios that SSD dominate their respective
sector indices. The proportion of the overall portfolio invested in each
sector is not pre-specified, rather it is decided via optimisation.

Computational results were given for our subset SSD approach as applied to the
S&P 500 over the period $29^{\text{th}}$ August 2018 to $29^{\text{th}}$
December 2023. These indicated that the scaled version of our subset SSD
approach significantly outperforms the S&P 500 over the period considered. Our
approach also outperforms the standard SSD based approach to the problem.

## References

  * Bawa [1975] V. S. Bawa.  Optimal rules for ordering uncertain prospects.  _Journal of Financial Economics_ , 2(1):95–121, 1975.  doi: 10.1016/0304-405X(75)90025-2.
  * Bruni et al. [2012] Renato Bruni, Francesco Cesarone, Andrea Scozzari, and Fabio Tardella.  A new stochastic dominance approach to enhanced index tracking problems.  _Economics Bulletin_ , 32(4):3460–3470, 2012.
  * Bruni et al. [2017] Renato Bruni, Francesco Cesarone, Andrea Scozzari, and Fabio Tardella.  On exact and approximate stochastic dominance strategies for portfolio selection.  _European Journal of Operational Research_ , 259(1):322–329, 2017.  doi: 10.1016/j.ejor.2016.10.006.
  * Cesarone and Puerto [2024] Francesco Cesarone and Justo Puerto.  New approximate stochastic dominance approaches for enhanced indexation models.  https://arxiv.org/html/2401.12669v19, 2024.
  * Cesarone et al. [2023] Francesco Cesarone, Raffaello Cesetti, Giuseppe Orlando, Manuel Luis Martino, and Jacopo Maria Ricci.  Comparing SSD-efficient portfolios with a skewed reference distribution.  _Mathematics_ , 11(1), 2023.  doi: 10.3390/math11010050.
  * CPLEX Optimizer 22.1.0 [2023] CPLEX Optimizer 22.1.0.  IBM. Available from https://www.ibm.com/products/ ilog-cplex-optimization-studio/cplex-optimizer/, last accessed October 17th 2023, 2023.
  * Dentcheva and Ruszczynski [2003] Darinka Dentcheva and Andrzej Ruszczynski.  Optimization with stochastic dominance constraints.  _SIAM Journal on Optimization_ , 14(2):548–566, 2003.  doi: 10.1137/S1052623402420528.
  * Dentcheva and Ruszczynski [2006] Darinka Dentcheva and Andrzej Ruszczynski.  Portfolio optimization with stochastic dominance constraints.  _Journal of Banking & Finance_, 30(2):433–451, 2006.  doi: 10.1016/j.jbankfin.2005.04.024.
  * Fábián et al. [2011a] C. Fábián, G. Mitra, and D. Roman.  Processing second-order stochastic dominance models using cutting-plane representations.  _Mathematical Programming_ , 130(1):33–57, 2011a.  doi: 10.1007/s10107-009-0326-1.
  * Fábián et al. [2011b] C. Fábián, G. Mitra, D. Roman, and V. Zverovich.  An enhanced model for portfolio choice with SSD criteria: a constructive approach.  _Quantitative Finance_ , 11(10):1525–1534, 2011b.  doi: 10.1080/14697680903493607.
  * Goel and Sharma [2021] Anubha Goel and Amita Sharma.  Deviation measure in second-order stochastic dominance with an application to enhanced indexing.  _International Transactions in Operational Research_ , 28(4):2218–2247, 2021.  doi: 10.1111/itor.12629.
  * Hadar and Russell [1969] J. Hadar and W. Russell.  Rules for ordering uncertain prospects.  _The American Economic Review_ , 59(1):25–34, 1969.  doi: 10.2307/1811090.
  * Hodder et al. [2015] James E Hodder, Jens Carsten Jackwerth, and Olga Kolokolova.  Improved portfolio choice using second-order stochastic dominance.  _Review of Finance_ , 19(4):1623–1647, 2015.  doi: 10.1093/rof/rfu025.
  * Kopa and Post [2015] M. Kopa and T. Post.  A general test for SSD portfolio efficiency.  _OR Spectrum_ , 37(1):703–734, 2015\.  doi: 10.1007/s00291-014-0373-8.
  * Künzi-Bay and Mayer [2006] Alexandra Künzi-Bay and János Mayer.  Computational aspects of minimizing conditional value-at-risk.  _Computational Management Science_ , 3(1):3–27, 2006.  doi: 10.1007/s10287-005-0042-0.
  * Kuosmanen [2001] Timo Kuosmanen.  Stochastic dominance efficiency tests under diversification.  https://econwpa.ub.uni-muenchen.de/econ-wp/fin/papers/0105/0105001.pdf, 2001\.
  * Kuosmanen [2004] Timo Kuosmanen.  Efficient diversification according to stochastic dominance criteria.  _Management Science_ , 50(10):1390–1406, 2004.  doi: 10.1287/mnsc.1040.0284.
  * Levy [1992] Haim Levy.  Stochastic dominance and expected utility: survey and analysis.  _Management Science_ , 38(4):555–593, 1992.  doi: 10.1287/mnsc.38.4.555.
  * Liesio et al. [2023] Juuso Liesio, Markku Kallio, and Nikolaos Argyris.  Incomplete risk-preference information in portfolio decision analysis.  _European Journal of Operational Research_ , 304(3):1084–1098, FEB 1 2023.  ISSN 0377-2217.  doi: 10.1016/j.ejor.2022.04.043.
  * Liu et al. [2021] Jia Liu, Zhiping Chen, and Giorgio Consigli.  Interval-based stochastic dominance: theoretical framework and application to portfolio choices.  _Annals of Operations Research_ , 307(1):329–361, 2021.  doi: 10.1007/s10479-021-04231-9.
  * Lizyayev and Ruszczyński [2012] Andrey Lizyayev and Andrzej Ruszczyński.  Tractable almost stochastic dominance.  _European Journal of Operational Research_ , 218(2):448–455, 2012.  doi: 10.1016/j.ejor.2011.11.019.
  * Luedtke [2008] James Luedtke.  New formulations for optimization under stochastic dominance constraints.  _SIAM Journal on Optimization_ , 19(3):1433–1450, 2008.  doi: 10.1137/070707956.
  * Malavasi et al. [2021] Matteo Malavasi, Sergio Ortobelli Lozza, and Stefan Truck.  Second order of stochastic dominance efficiency vs mean variance efficiency.  _European Journal of Operational Research_ , 290(3):1192–1206, MAY 1 2021.  ISSN 0377-2217.  doi: 10.1016/j.ejor.2020.08.051.
  * Markowitz [1952] H. Markowitz.  Portfolio selection.  _Journal of Finance_ , 7(1):77–91, 1952\.  doi: 10.1111/j.1540-6261.1952.tb01525.x.
  * Ogryczak and Ruszczynski [2002] W. Ogryczak and A. Ruszczynski.  Dual stochastic dominance and related mean-risk models.  _SIAM Journal on Optimization_ , 13(1):60–78, 2002.  doi: 10.1137/S1052623400375075.
  * Post and Kopa [2013] T. Post and M. Kopa.  General linear formulations of stochastic dominance criteria.  _European Journal of Operational Research_ , 230(2):321–332, 2013.  doi: 10.1016/j.ejor.2013.04.015.
  * Post [2003] Thierry Post.  Empirical tests for stochastic dominance efficiency.  _Journal of Finance_ , 58(5):1905–1931, 2003.  doi: 10.1111/1540-6261.00592.
  * Roman et al. [2006] D. Roman, K. Darby-Dowman, and G. Mitra.  Portfolio construction based on stochastic dominance and target return distributions.  _Mathematical Programming_ , 108(2):541–569, 2006.  doi: 10.1007/s10107-006-0722-8.
  * Roman et al. [2013] D. Roman, G. Mitra, and V. Zverovich.  Enhanced indexation based on second-order stochastic dominance.  _European Journal of Operational Research_ , 228(1):273–281, 2013.  doi: 10.1016/j.ejor.2013.01.035.
  * Sehgal and Mehra [2021] Ruchika Sehgal and Aparna Mehra.  Robust reward–risk ratio portfolio optimization.  _International Transactions in Operational Research_ , 28(4):2169–2190, 2021.  doi: 10.1111/itor.12652.
  * Sharma and Mehra [2017] Amita Sharma and Aparna Mehra.  Financial analysis based sectoral portfolio optimization under second order stochastic dominance.  _Annals of Operations Research_ , 256:171–197, 2017\.  doi: 10.1007/s10479-015-2095-y.
  * Sharma et al. [2017] Amita Sharma, Shubhada Agrawal, and Aparna Mehra.  Enhanced indexing for risk averse investors using relaxed second order stochastic dominance.  _Optimization and Engineering_ , 18(2):407–442, 2017.  doi: 10.1007/s11081-016-9329-y.
  * Valle et al. [2017] C. A. Valle, D. Roman, and G. Mitra.  Novel approaches for portfolio construction using second order stochastic dominance.  _Computational Management Science_ , 14(2):257–280, 2017.  doi: 10.1007/s10287-017-0274-9.
  * Whitmore and Findlay [1978] G. A. Whitmore and M. C. Findlay.  _Stochastic dominance: an approach to decision-making under risk_.  Lexington Books, 1978.


# Decentralized Multi-Agent Planning for Multirotors:
a Fully Online and Communication Latency Robust Approach

Charbel Toumieh The author is an independent researcher (e-mail:
<EMAIL_ADDRESS>
###### Abstract

There are many industrial, commercial and social applications for multi-agent
planning for multirotors such as autonomous agriculture, infrastructure
inspection and search and rescue. Thus, improving on the state-of-the-art of
multi-agent planning to make it a viable real-world solution is of great
benefit. In this work, we propose a new method for multi-agent planning in a
static environment that improves our previous work by making it fully online
as well as robust to communication latency. The proposed framework generates a
global path and a Safe Corridor to avoid static obstacles in an online fashion
(generated offline in our previous work). It then generates a time-aware Safe
Corridor which takes into account the future positions of other agents to
avoid intra-agent collisions. The time-aware Safe Corridor is given with a
local reference trajectory to an MIQP (Mixed-Integer Quadratic Problem)/MPC
(Model Predictive Control) solver that outputs a safe and optimal trajectory.
The planning frequency is adapted to account for communication delays. The
proposed method is fully online, real-time, decentralized, and synchronous. It
is compared to 3 recent state-of-the-art methods in simulations. It
outperforms all methods in robustness and safety as well as flight time. It
also outperforms the only other state-of-the-art latency robust method in
computation time.

video: https://youtu.be/eKwYNU1Q0wY

##  I INTRODUCTION

###  I-A Problem statement

Multi-agent planning has been gaining in popularity in the research community
due to recent advances. These advances are making it a viable solution to many
commercial, industrial, and military applications. There are multiple
challenges that face a multi-agent planning framework such as the problem of
synchronizing agents for synchronous planning methods and dealing with
communication latency. It is the purpose of this paper to extend upon our
previous state-of-the-art work [1] that outperformed other state-of-the-art
methods in computation efficiency, trajectory speed, and smoothness in a
cluttered environment. We provide a new approach derived from [1] that is
fully online and robust to arbitrary communication latency. We also study the
effect of communication latency on the overall performance of our planner and
compare it with other state-of-the-art methods.

###  I-B Related work

####  I-B1 Multi-agent planning for multirotors

In [2], the authors present a centralized multi-agent planning framework that
uses time-aware Safe Corridors. The method has 3 sequential steps: roadmap
generation, then discrete planning, and finally continuous refinement. The
approach presented by the authors is centralized although some steps can be
decentralized. While the computation time is not suitable for online high-
speed planning and replanning, the method used served as an inspiration for
many subsequent methods in the state-of-the-art. Such methods include [3] and
[1] which in turn served as an inspiration for the work presented in this
paper.

Buffered Voronoi Cells have been used by multiple works [4], [5] for multi-
agent collision avoidance but do not account for static obstacles. Other
approaches [6] use separating hyperplanes to avoid collisions between agents
and model static obstacles in the form of ellipsoid constraints in a
decentralized MPC formulation. The generation of ellipsoid representation of
the environment is not trivial and is not addressed by the authors of [6].

MADER, an asynchronous multi-agent planning framework has been proposed in
[7]. The method allows for avoiding static, and dynamic obstacles, as well as
other planning agents. The authors combine a search-based approach with an
optimization approach, where the output of the search-based approach is taken
as initialization for the optimization problem. This choice was made since the
optimization problem defined by the authors is non-convex and requires a good
initial guess.

EGO-Swarm was proposed in [8] as an asynchronous and decentralized trajectory
planner. It requires each planning agent to broadcast its generated trajectory
at a fixed frequency. When each agent receives the trajectories of other
agents, it proceeds immediately to do a collision check. While the approach
has been demonstrated in real-world experiments, it still suffers from
collisions due to communication delays between agents.

In a similar fashion to [2], the authors of [3] present a distributed and
online trajectory generation framework for multi-agent quadrotor systems using
time-aware Safe Corridors (or Linear Safe Corridors). The environment
representation used by the authors is an octomap [9]. The Safe Corridor used
to generate the time-aware Safe Corridor contains only one polyhedron which
leads to slow and conservative trajectories.

In [10], a decentralized model predictive control approach is used for
collision avoidance and cohesive flight. The obstacles are described as
mathematical functions (cylinders, paraboloids …) in order to include them in
the decentralized MPC formulation as constraints. It is however not trivial to
describe an arbitrary cluttered environment through continuous mathematical
functions that are easy to add as constraints to an MPC formulation.

Finally, in our previous work [1], we proposed a decentralized and synchronous
planning framework that is inspired by [2]. The approach takes into account
static obstacles using Safe Corridors (generated from a voxel grid
representation [11]). Safe Corridors are then augmented to time-aware Safe
Corridors to avoid intra-agent collisions. The proposed approach outperforms
state-of-the-art methods in all performance metrics, including robustness,
computation time, and trajectory speed.

####  I-B2 Latency robust multi-agent planning

The previously cited works do not account for communication delay, or can
passively handle latency up to a fixed limit [1]. Some multi-agent planning
frameworks take into account communication delay and will be presented in this
section.

In [12], an asynchronous and decentralized trajectory planner is presented.
The planner guarantees safety using separating hyperplanes from previous
planning iterations. While the presented approach can handle communication
delays, it does not account for any type of obstacles (static or dynamic),
which limits its applicability to the real world.

Finally, RMADER (Robust MADER) is proposed in [13], which is an extension of
MADER [7]. They convexify the optimization problem in order to improve the
computation time. However, they inherit from MADER the polyhedral
representation of the obstacles in the environment. This representation is not
trivial to generate and can add significant overhead to the planning
framework.

###  I-C Contribution

The main contribution of our paper is an improved decentralized and
synchronous planning framework that is robust to communication latency. The
proposed framework is built on our previous work [1] and conserves its
advantages. Thus, the proposed method has low computation time and takes into
account static obstacles and other planning agents. The improvements are:

  1. 1.

The addition of a mechanism to deal with arbitrary communication latency by
dynamically adapting the planning frequency to avoid collisions and guarantee
safety.

  2. 2.

The integration of 2 previously offline steps in [1] (global path generation
step and Safe Corridor generation step) to make the framework fully online and
suitable for real-world applications.

  3. 3.

The modification of the stalemate/deadlock resolution mechanism to guarantee
safety.

The method is tested in simulations to show the effect of communication
latency on the performance of the planner. It is also compared to 3 recent
works: EGO-Swarm [8], MADER [7] and RMADER [13] in terms of trajectory
safety/performance as well as computation time.

##  II Assumptions

Figure 1: We show the global pipeline of the planning framework of a single
planning agent. It is run in a loop at a varying/adaptive frequency.

We assume perfect control (the controller executes the generated trajectory
perfectly) and perfect localization (each agent can localize itself and other
agents at any moment to an arbitrary accuracy). These assumptions are made by
all of the previously cited state-of-the-art methods. In addition to these
assumptions, we assume that the clocks of the agents are synchronized. We
assume 2 cases:

  1. 1.

We can synchronize all agents at the beginning of a given mission.

  2. 2.

If an agent (not synchronized) is getting close to a cluster of other
synchronized agents, we assume the range of communication is big enough so
that the agent can synchronize its clock with the cluster before getting close
enough for collision avoidance.

Furthermore, we assume symmetric behavior of the communication: if there is a
latency in the delivery of a message from agent $i$ to agent $j$ in a given
planning iteration/period, the same latency happens when agent $j$ is trying
to deliver a message to agent $i$.

(a) Safe Corridor at iteration $k$.

(b) Safe Corridor at iteration $k+1$.

Figure 2: The obstacles are shown in red. The predicted positions of the agent
are shown as yellow circles (MPC trajectory). They get increasingly
transparent as we move forward in time. At iteration $k$ (Fig. 2(a)), all
polyhedra (in blue) contain at least one point of the MPC trajectory. At the
next iteration $k+1$ (Fig. 2(b)), the first position of the MPC trajectory
moves out of the first polyhedron (in dashed blue lines). Thus, we remove it
from the Safe Corridor and generate another polyhedron (in green) using the
global path. The new polyhedron is added to the Safe Corridor.

##  III The planner

Our planner is run concurrently on each agent in a swarm. The dynamical model
of each agent is the same as presented in [1]. We use a voxel grid
representation of the environment, which can be trivially and efficiently
generated [11]. Each agent has a voxel grid that is of fixed size and that
moves with the agent such that the agent is always at its center. This voxel
grid is used for global path finding and Safe Corridor generation. The clocks
of the agents are synchronized.

In [1], the planning is divided into 2 stages: an offline stage for global
path finding and Safe Corridor generation; then an online stage where the
time-aware Safe Corridors and the dynamically feasible trajectory are
generated. In the planner proposed in this paper, the offline stage is now
integrated into the online planning stage so the whole planning/replanning
framework is run online. This makes it suitable for real-world deployment and
missions such as exploration. The steps of the proposed planner are (Fig. 1):

  1. 1.

Generate a global path (Sect. III-A).

  2. 2.

Generate a Safe Corridor (Sect. III-B)

  3. 3.

Generate a time-aware Safe Corridor (Sect. III-C).

  4. 4.

Generate a local reference trajectory (Sect. III-D).

  5. 5.

Solve the Mixed-Integer Quadratic Program (MIQP)/Model Predictive Control
(MPC) problem to generate a locally optimal trajectory (Sect. III-E).

In the first step, we generate a global path from the position of the agent to
the goal position. This path avoids all static obstacles and is used to
generate the Safe Corridor and to generate the local reference trajectory. In
the second step, we generate a Safe Corridor (a series of overlapping convex
polyhedra) that covers only the free space in the environment. These convex
polyhedra are used as linear constraints in an optimization formulation to
constrain the trajectory to the free space and avoid collisions with static
obstacles. In the third step, we use the recently generated trajectories of
the agents and the Safe Corridor to generate time-aware Safe Corridors. This
allows the agents to avoid intra-agent collisions. In the fourth step, we
sample the global path at a given velocity to generate a local reference
trajectory that the dynamically feasible trajectory tries to follow as closely
as possible. In the fifth and final step, we generate the dynamically feasible
trajectory to be executed by the agent. It is generated by solving an
optimization problem that takes time-aware Safe Corridors and a local
reference trajectory and guarantees that there are no collisions of any nature
(intra-agent or static obstacles) while the agent moves closer to its goal.

These steps were run sequentially and periodically at a fixed frequency in our
previous work [1]. However, in this work, we vary the planning frequency to
account for communication latency. As in [1], each agent broadcasts its
planned trajectory at the end of the planning iteration so that other agents
can know it. In addition to the planned trajectory, we also broadcast the
times we started and finished generating the trajectory so that other agents
can estimate the communication latency (not done in [1] \- more details in
Sect. III-F). We briefly explain each step in this section while focusing more
on the steps where changes were made with respect to [1].

###  III-A Generate a global path

In this step, a global path is generated connecting the current position of
the agent to the desired final position using the local voxel grid. The
occupied voxels in the voxel grid are inflated by each agent’s size before
feeding the grid to the path planning algorithm. In case the goal position is
outside the local voxel grid of the agent, we choose an intermediate goal in
the grid as presented in [14]. The main idea is to draw a line connecting the
position of the agent to the goal and get the intersection with the borders of
the voxel grid. This intersection is a voxel and is set as an intermediate
goal. We also clear/set to free all the border voxels of the voxel grid to
help the agent find a path to the intermediate goal in extremely cluttered
environments.

At each iteration, the starting point for the global path search is the last
point in the local reference trajectory generated in the previous planning
iteration (Sect. III-D). The local reference trajectory is then connected to
the path found through the global search to generate the final global path
used in the subsequent sections (for generating the local reference trajectory
of the current iteration).

We use JPS (Jump Point Search) [15] and DMP (Distance Map Planner) for path
planning. JPS employs pruning techniques on the A* algorithm to potentially
speed up the generation time by an order of magnitude. DMP uses artificial
potential fields to push the path generated by JPS away from obstacles. This
adds an additional margin of safety and improves the trajectory generated in
the last step (MIQP optimization output) in terms of speed and smoothness (see
[14] for more details).

(a) Stalemate caused by a symmetrical position.

(b) Perturbing hyperplanes asymmetrically.

(c) Perturbing hyperplanes symmetrically.

Figure 3: A stalemate/deadlock happens when 2 agents are trying to move
towards opposite goals and the solver is stuck on the borders of the
hyperplanes (Fig. 3(a)). Any movement up or down would not decrease the
distance to the goal. If the hyperplanes are perturbed asymmetrically as done
in [1] (Fig. 3(b)), the distance between the agents can potentially become
lower than the safety distance. We modify the perturbation vector (Sect.
III-C) to make the perturbation symmetrical and guarantee safety when the
agents move in the direction of the magenta vectors or any other direction
(Fig. 3(c)). Figure 4: We show the trajectories of 2 agents (in red and
yellow) and the corresponding discrete positions that get more transparent as
we move forward in time. We ignore the positions of each trajectory that have
no corresponding position in the other ($k-2$ and $k+3$). The separating
hyperplanes (dashed lines in different colors) are generated between the
positions of the agents corresponding to the same time in the future starting
from the current iteration $k$. The last separating hyperplane $k+2$ is used
to fill the remaining $N-3$ hyperplanes required to generate the TASC.

###  III-B Generate a Safe Corridor around the global path

Safe Corridors are a series of overlapping convex shapes that cover only free
space in the environment. They are used by many state-of-the-art planning
methods to constrain a dynamically feasible trajectory inside them, and thus
guarantee safety [16], [14], [1]. Many methods exist in the literature for
Safe Corridor generation [17], [18], [19] [20]. The method used for the
generation is [19] since it provides the best performance among the state-of-
the-art methods for trajectory planning.

The Safe Corridor generation method takes as input a voxel grid (the local
voxel grid centered around the agent) and the global path around which we want
to generate the Safe Corridor. At each iteration, we always make sure that we
have a certain number $P_{\text{hor}}$ of polyhedra that cover the free space
of the environment.

At the first iteration of planning, we use the global path at the first
iteration to generate a Safe Corridor that contains up to $P_{\text{hor}}$
number of polyhedra (polyhedra horizon). Subsequently, at each planning
period, we use the global path generated in this planning period to update the
Safe Corridor generated in the last step. The update consists of the following
(Fig. 2): all the polyhedra that contain at least one point of the last
generated MPC trajectory are kept. The other polyhedra are removed and new
polyhedra are generated in their place until we have $P_{\text{hor}}$
polyhedra in total. To generate each polyhedron, we sample the global path at
a constant step (voxel size). We then use the first point of the sampled
global path that is outside all the remaining polyhedra as a seed voxel to
generate an additional polyhedron.

###  III-C Generate a time-aware Safe Corridor (TASC)

After generating the Safe Corridor, we use it along with the trajectories
generated by all the other agents at the previous iterations to create a time-
aware Safe Corridor (TASC). The future positions predicted by the MPC
trajectories of the agents at the previous planning iterations are used to
generate hyperplanes to constrain the future/MPC positions at the current
iteration. These hyperplanes are added to the constraints of the Safe
Corridor. This creates a series of Safe Corridors at each planning iteration
that we call time-aware Safe Corridors in [1]. We refer the reader to [1] for
a detailed explanation of how time-aware Safe Corridors are generated.

We augment/improve the TASC generation method to account for trajectories that
were not generated at the same planning iteration $k$ (Fig. 4). We ignore the
positions of each trajectory that have no corresponding positions in the other
trajectory ($k-1$ and $k+3$ in Fig. 4). Then, starting with the position of
the current iteration $k$, we generate separating hyperplanes for the rest of
the common positions ($k$, $k+1$ and $k+2$ in Fig. 4). Since we need $N$
separating hyperplanes to generate the TASC (as shown in [1]), we set the rest
of the hyperplanes equal to the last separating hyperplanes ($k+2$ in Fig. 4).

####  III-C1 Dealing with stalemates/deadlocks

In [1], in order to avoid stalemates/deadlocks, we modified the normal vectors
of the separating hyperplanes by perturbing them constantly through time (a
time-varying right-hand rule). This would avoid adding an explicit mechanism
that creates subgoals for each agent to avoid stalemates/deadlocks like in
[21]. We defined the normalized plane normal
$\boldsymbol{n}_{\text{hyp,norm}}$, the right vector $\boldsymbol{r}$ that is
the cross product between $\boldsymbol{n}_{\text{hyp,norm}}$ and
$\boldsymbol{z}_{W}$ plus the cross product between
$\boldsymbol{n}_{\text{hyp,norm}}$ and $\boldsymbol{y}_{W}$, a perturbation
$m$, and a user-chosen coefficient $c$ that defines how tilted the final
normal vector of the hyperplane $\boldsymbol{n}_{\text{hyp,final}}$ is with
respect to the initial vector $\boldsymbol{n}_{\text{hyp}}$:

$\displaystyle\boldsymbol{n}_{\text{hyp,norm}}=\dfrac{\boldsymbol{n}_{\text{hyp}}}{||\boldsymbol{n}_{\text{hyp}}||_{2}}$
(1)
$\displaystyle\boldsymbol{z}_{W}=[0,0,1]^{T},\quad\boldsymbol{y}_{W}=[0,1,0]^{T}$
(2)
$\displaystyle\boldsymbol{r}=\boldsymbol{n}_{\text{hyp,norm}}\times\boldsymbol{z}_{W}+\boldsymbol{n}_{\text{hyp,norm}}\times\boldsymbol{y}_{W}$
(3)
$\displaystyle\boldsymbol{n}_{\text{pert}}=(c+m)\cdot\dfrac{\boldsymbol{r}}{||\boldsymbol{r}||_{2}}+c\cdot\boldsymbol{z}_{W}$
(4)
$\displaystyle\boldsymbol{n}_{\text{hyp,final}}=\boldsymbol{n}_{\text{pert}}+\boldsymbol{n}_{\text{hyp,norm}}$
(5)

However, a component of the perturbation vector $\boldsymbol{n}_{\text{pert}}$
is non-symmetric ($c\cdot\boldsymbol{z}_{W}$), which can generate normal
vectors that are non-colinear. This can result in cases where the distance
between agents is lower than the safety/collision distance $2\cdot
d_{\text{rad}}$ (Fig. 3). For this reason, we replace the non-symmetric term
with the following symmetric term:
$c\cdot(\boldsymbol{z}_{W}\times\boldsymbol{n}_{\text{hyp,norm}})$. The final
perturbation vector then becomes:

$\displaystyle\boldsymbol{n}_{\text{pert}}=(c+m)\cdot\dfrac{\boldsymbol{r}}{||\boldsymbol{r}||_{2}}+c\cdot(\boldsymbol{z}_{W}\times\boldsymbol{n}_{\text{hyp,norm}})$
(6)

It is then added to $\boldsymbol{n}_{\text{hyp,norm}}$ to generate
$\boldsymbol{n}_{\text{hyp,final}}$ as in equation (5).

###  III-D Generate a local reference trajectory

We use the global path to generate a local reference trajectory that is used
as a reference for the MPC to follow. The generation of such reference
trajectory is done by sampling the global path at a constant velocity
$v_{\text{samp}}$. The number of sampled points is equal to the number of
discretization steps ($N$) in the MPC/MIQP formulation.

We only generate a new local reference trajectory in the following case: the
last point of the MPC trajectory is within a distance $d_{\text{thresh}}$ from
the last point of the local reference trajectory generated at the previous
iteration. Otherwise, we keep the local reference trajectory generated at the
previous planning iteration.

Figure 5: We show an example of how different agents handle communication
delays between each other. In this example agent 2 communicates with agents 1
and 3, whereas agents 1 and 3 do not communicate with each other (not within
the range of communication). We show in green the computation time of each
agent, in blue the communication latency between agents 1 and 2, and in red
the communication latency between agents 2 and 3. The arrows indicate the time
at which an agent $i$ receives the trajectory $\boldsymbol{T}_{j,k}$ of
another agent $j$ generated at iteration $k$. At the first iteration, all
agents synchronize their first planning iteration to be at the same time. At
the subsequent iterations, an agent skips planning in one of 2 cases: 1) At
least one agent within the communication range is yet to receive its last
generated trajectory 2) It is yet to receive a new generated trajectory of
another agent within the communication range and it has used all the
previously received trajectories of this agent to generate its own trajectory.

###  III-E Solving the MIQP/MPC problem

In this final step, we take the reference trajectory, and we solve an MPC
optimization problem that minimizes the distance of the generated trajectory
to the reference trajectory while also minimizing the jerk for smoothness. The
generated trajectory consists of $N+1$ discrete states $\boldsymbol{x}_{i}$,
$i=0,1,...,N$ that contain the position, velocity, and acceleration of the
agent. Each consecutive pair of discrete states are separated by a time step
$h$. Thus, the time horizon of the planning is $N\cdot h$. The velocity and
acceleration of the last state $\boldsymbol{x}_{N}$ are constrained/set to 0
to guarantee a safe trajectory for all agents in case subsequent optimizations
fail (see [1] for more details).

The time-aware Safe Corridor is used to ensure the safety of the trajectory.
We add the linear constraints of the time-aware Safe Corridor to the MPC
optimization problem. By forcing each segment of the MPC trajectory be in at
least one of the polyhedra of the time-aware Safe Corridor, we ensure no
collision happens between the agent and the static obstacles as well as other
planning agents. The final formulation of the optimization problem is a Mixed-
Integer Quadratic Problem (MIQP) exactly like the one presented in [14], [1].

###  III-F Handling communication delay

Algorithm 1 Run at every iteration $k$ for agent $i$:

1:delay_planning = false

2:for each agent $j$ in $J$ do

3: if received $\boldsymbol{T}_{j}$ then

4: traj_old[$j$].add($\boldsymbol{T}_{j}$)

5: else

6: if traj_old[$j$].size() == 0 then

7: delay_planning = true

8: if not(delay_planning) then

9: $dt_{\text{delay},i,j}$ = ComputeLatency(traj_old[$j$][0])

10: if $dt_{\text{delay},i,j}+\boldsymbol{T}_{i,\text{last}}$.end then
$>t_{\text{cur}}$

11: delay_planning = true

12:if not(delay_planning) then

13: for each agent $j$ in $J$ do

14: GenerateTASC(traj_old[$j$][0], $\boldsymbol{T}_{i,\text{last}}$)

15: traj_old[$j$].RemoveFirstElement()

Our previous work [1] ran the planning algorithm at a constant period equal to
the MPC discretization step $dt_{\text{plan}}=h$. It was able to handle
communication delay passively by assuming that the communication delay was
lower than a time variable $dt_{\text{max,delay}}$ equal to the planning
period $dt_{\text{plan}}$ minus the planner computation time
$dt_{\text{comp}}$
($dt_{\text{max,delay}}=dt_{\text{plan}}-dt_{\text{comp}}$). However, no
mechanism was in place to handle the communication latency when it exceeds
$dt_{\text{max,delay}}$.

In this work, we propose to adapt the planning period to be able to guarantee
safety no matter the communication delay. In addition to broadcasting the
trajectory $\boldsymbol{T}_{j}$ when it finishes generating it, each agent $j$
broadcasts the time at which it started generating its trajectory i.e. the
time at the start of the planning period ($\boldsymbol{T}_{j}$.start). It also
broadcasts the time it finished generating the trajectory i.e. the time it
sent it ($\boldsymbol{T}_{j}$.end). This allows another agent $i$ to estimate
the communication delay between it and agent $j$ since their clocks are
synchronized. The delay can be estimated by subtracting
$\boldsymbol{T}_{j}$.end from the reception time of agent $i$,
$t_{\text{rec},i}$:

$\displaystyle
dt_{\text{delay},i,j}=t_{\text{rec},i}-\boldsymbol{T}_{j}\text{.end}$ (7)

This in turn allows agent $i$ to know whether its last generated trajectory
$\boldsymbol{T}_{i,\text{last}}$ was received by agent $j$ before the start
time of the current planning period $t_{\text{cur}}$. The last generated
trajectory of agent $i$ is not yet received by agent $j$ if the following
condition is true:

$\displaystyle
dt_{\text{delay},i,j}+\boldsymbol{T}_{i,\text{last}}.\text{end}>t_{\text{cur}}$
(8)

The planner will skip planning at the start of the current planning period and
wait for the next period if one of these 2 cases is true:

  1. 1.

It knows that there is another agent within its communication range that is
yet to receive its last planned trajectory.

  2. 2.

It is yet to receive a new planned trajectory of another agent within its
communication range and it has used all the old received trajectories of this
agent for planning.

We propose the following algorithm to handle communication latency (Alg. 1).
At every planning iteration (which happens every $dt_{\text{plan}}=h$), every
agent $i$ checks if it received a trajectory from every other agent $j$ (line
3). If it did, it adds the received trajectory to a 2D vector (traj_old) whose
first index indicates the number or ID of the other agent i.e. $j$ (line 4).
If agent $i$ did not receive a trajectory from agent $j$, it checks if there
is an unused old trajectory in the vector traj_old[$j$] (line 5-6). If not, we
delay the planning since we have no new or old trajectory to use for
generating the TASC (line 7). If the planning should not be delayed due to
previous conditions (line 8), we check if it should be delayed because agent
$j$ hasn’t received the trajectory of agent $i$ yet. This is done by first
computing the communication delay using equation (7) (line 9), and then
checking the condition (8) (lines 10-11). Finally, we check if the planning
should be delayed after going through all agents (line 12). If not, we compute
the TASC using the oldest unused trajectory of each agent $j$ and remove it
from the vector of old trajectories (lines 13-15). The starting time
$\boldsymbol{T}_{j}$.start allows to know at which iteration $k$ the
trajectory was generated, which is important in TASC generation (Fig. 4).

We show an example of how this algorithm would perform in Fig. 5. In this
example, agent 2 sees and communicates with agents 1 and 3, but agents 1 and 3
do not see and communicate with each other. Still, the algorithm allows for
safe planning and coordination between all agents.

##  IV Simulation Results

The testing setup is similar to what is presented in [13]. Thus, we will use
their results as a reference for our comparison. The simulations are run on
Intel i7 CPUs with a base frequency of 2.6GHz and a turbo boost of 4GHz. The
testing consists of 10 agents in a circular configuration (Fig. 6(a))
exchanging positions. We compare our method with RMADER [13] and 2 versions of
Ego-Swarm [8]. We set the maximum velocity $v_{\text{max}}=10\ \text{m/s}$,
the maximum acceleration $a_{\text{max}}=20\ \text{m/s\textsuperscript{2}}$
and the maximum jerk $j_{\text{max}}=30\ \text{m/s\textsuperscript{3}}$ for
RMADER, Ego-Swarm and our method (along the $x$, $y$ and $z$ directions). For
Ego-Swarm, we also consider a more conservative version (slow Ego-Swarm) with
a maximum acceleration $a_{\text{max}}=10\ \text{m/s\textsuperscript{2}}$ and
a maximum velocity $v_{\text{max}}=5\ \text{m/s}$.

For MADER and RMADER, each agent is represented as a bounding box of size
$0.25\times 0.25\times 0.25$ m. For Ego-Swarm and our planner, each agent is
represented as a sphere of diameter $0.25$ m as per the experiments in [13]
(at the time of writing, the bounding box dimensions and sphere diameter were
not mentioned in [13], but they were communicated to us by the authors of
[13]). The comparison is done with 100 simulated runs for communication
latencies equal to $0$, $50$, and $100$ milliseconds. The comparison metrics
are:

  1. 1.

Collision %: percentage of simulations where there was at least one collision.

  2. 2.

Average number of stops expected in a single simulation from all agents.

  3. 3.

Mean of the jerk cost
$J_{\text{cost}}=\int_{t_{\text{ini}}}^{t_{\text{fin}}}||\boldsymbol{j}(t)||^{2}\mathrm{d}t$
where $t_{\text{ini}}$ and $t_{\text{fin}}$ are the initial and final time of
the trajectory.

  4. 4.

Mean of the acceleration cost
$A_{\text{cost}}=\int_{t_{\text{ini}}}^{t_{\text{fin}}}||\boldsymbol{a}(t)||^{2}\mathrm{d}t$.

  5. 5.

Mean and max flight time.

  6. 6.

Computation time.

Table I: Comparison between Ego-Swarm (ES) [8], slow Ego-Swarm (Slow ES) [8], MADER [7], RMADER [13] and our method. The comparison consists of 100 simulations with communication delays between 10 agents exchanging positions in a circular configuration as in Fig. 6(a). The communication delays are $dt=0$ ms $\mid$ $dt=50$ ms $\mid$ $dt=100$ ms. We show in bold the best performer among the safe planners (RMADER [13] and our planner). Method | Collision [%] | Mean # stops | Accel. cost (m/s2) | Jerk cost (103 m/s3) | Mean flight time (s) | Max flight time (s)
---|---|---|---|---|---|---
ES [8] |  64 $\mid$ 84 $\mid$ 84 |  0.004 $\mid$ 0 $\mid$ 0.01 |  662 $\mid$ 700 $\mid$ 788 |  9.07 $\mid$ 9.46 $\mid$ 10.4 |  7.19 $\mid$ 7.24 $\mid$ 7.28 |  7.38 $\mid$ 7.51 $\mid$ 7.63
Slow ES [8] |  14 $\mid$ 25 $\mid$ 22 |  0 $\mid$ 0 $\mid$ 0 |  110 $\mid$ 113 $\mid$ 113 |  15.4 $\mid$ 15.5 $\mid$ 15.5 |  11.6 $\mid$ 11.7 $\mid$ 11.8 |  11.9 $\mid$ 12 $\mid$ 13
MADER [7] |  15 $\mid$ 38 $\mid$ 42 |  0 $\mid$ 0.001 $\mid$ 0 |  78.1 $\mid$ 74.2 $\mid$ 74.5 |  1.59 $\mid$ 1.64 $\mid$ 1.64 |  6.28 $\mid$ 6.25 $\mid$ 6.26 |  7.15 $\mid$ 7.35 $\mid$ 7.04
RMADER [13] |  0 $\mid$ 0 $\mid$ 0 |  0.46 $\mid$ 0.347 $\mid$ 1.75 |  127 $\mid$ 148 $\mid$ 190 |  2.94 $\mid$ 3.71 $\mid$ 5.94 |  7.28 $\mid$ 7.95 $\mid$ 10.4 |  8.41 $\mid$ 8.80 $\mid$ 11.9
proposed |  0 $\mid$ 0 $\mid$ 0 |  0 $\mid$ 0 $\mid$ 0 |  109 $\mid$ 114 $\mid$ 119 |  2.27 $\mid$ 2.49 $\mid$ 5.03 |  6.77 $\mid$ 6.79 $\mid$ 7.1 |  7.1 $\mid$ 7.3 $\mid$ 7.7

(a) Our planner: 10 agents with $dt=100$ ms with the setup in Tab. I.

(b) Our planner: 12 agents with $dt=0$ ms and obstacles (Sect. IV-C).

(c) Our planner: 12 agents with $dt=150$ ms and obstacles (Sect. IV-C).

Figure 6: The agents start in a circular configuration and swap positions. We show an overhead view of the trajectories generated by our planner in different settings (with and without obstacles), different communication latencies, and different dynamic limits. Table II: Computation time of our planner for the results in Tab. I. We show the mean / max / standard deviation. |  $dt=0$ ms |  $dt=50$ ms |  $dt=100$ ms
---|---|---|---
Comp. (ms) | 10.4 / 61 / 6.6 | 10.1 / 54.7 / 6.4 | 11.4 / 70 / 6.7

###  IV-A Planner parameters

The local voxel grid around each agent is of size $15\times 15\times 3.3$ m
and has a voxel size of $0.3$ m. We choose the following parameters: $N=9$,
$h=100$ ms, $v_{\text{samp}}=4.5$ m/s, $P_{\text{hor}}=3$,
$d_{\text{thresh}}=0.4$ m. The rest of the parameters are chosen the same as
in [1] with the exception of the maximum velocity, acceleration, and jerk
which are the same for all planners (Sect. IV).

###  IV-B Comparison with the state-of-the-art

We show in Tab. I the results of the planners with different communication
latencies ($0$, $50$, and $100$ ms). Our planner and Ego-Swarm [8] use voxel
girds as representations of the obstacles in the environment. MADER [7] and
RMADER [13] on the other hand use a polyhedral representation of the
environment i.e. all obstacles are represented by a series of convex
polyhedra. This representation is not trivial to generate and may add
considerable overhead to the autonomous navigation pipeline.

Our planner and RMADER [13] are the only planners that are able to generate
collision-free trajectories in all simulations, so we will focus our
comparison on them. Our planner outperforms RMADER in trajectory smoothness
across all latencies using both the acceleration ($25$% better on average) and
the jerk ($24$% better on average) metrics.

The mean and max flight times of our planner grow slower than those of RMADER
with the increase in latency. Over all latencies, our planner outperforms
RMADER in mean flight time by an average of $18$% and max flight time by an
average of $23$%.

####  IV-B1 Computation time

Ego-Swarm is the most computationally efficient with an average computation
time of $0.5$ ms. RMADER improves on MADER [7] in computation time by changing
the optimization problem from non-convex to convex. This improves the mean
computation time by $20$% (from $39.23$ ms to $31.08$ ms) and the max
computation time by $40$% (from $724$ ms to $433$ ms) as reported in [13].
While our planner is not as efficient as Ego-Swarm, it is much more efficient
than RMADER as shown in Tab. II. The mean computation time across all
latencies is $10.6$ ms and the max is $70$ ms.

Table III: Results for 8 and 12 agents in an environment with obstacles (Sect. IV-C). The mean / max / standard deviation of each metric is shown. # |  $dt$ (ms) | Distance (m) | Velocity (m/s) | Flight time (s) | Comp. time (ms) | Acc. cost (m/s2) | Jerk cost (103m/s3)
---|---|---|---|---|---|---|---
8 | 0 | 21.6 / 23.1 / 0.72 | 2.52 / 4.21 / 1.24 | 8.47 / 9.5 / 0.4 | 5.5 / 48.7 / 3 | 121 / 170 / 26.2 | 3.5 / 5.56 / 0.94
50 | 21.6 / 23.1 / 0.72 | 2.51 / 4.21 / 1.24 | 8.47 / 9.5 / 0.4 | 5.4 / 48.1 / 2.9 | 121 / 170 / 26.2 | 3.5 / 5.56 / 0.95
100 | 21.6/ 23.4 / 0.76 | 2.43 / 4.24 / 1.22 | 8.7 / 9.5 / 0.42 | 6.2 / 35 / 3.8 | 124 / 182 / 26.7 | 6.59 / 9.11 / 0.96
150 | 21.6 / 23.4 / 0.76 | 2.43 / 4.24 / 1.22 | 8.7 / 9.5 / 0.42 | 6.1 / 33.3 / 3.8 | 124 / 182 / 26.5 | 6.59 / 9.11 / 0.96
12 | 0 | 21.7 / 24.2 / 0.73 | 2.45 / 4.5 / 1.23 | 8.7 / 9.9 / 0.45 | 8.7 / 72.4 / 6 | 130 / 207 / 26.8 | 3.76 / 5.65 / 0.84
50 | 21.7 / 24.1 / 0.73 | 2.46 / 4.5 / 1.24 | 8.7 / 9.9 / 0.43 | 8.4 / 69.6 / 5.8 | 136 / 207 / 27.6 | 4.19 / 6.56 / 0.88
100 | 21.6 / 23.9 / 0.71 | 2.38 / 4.36 / 1.2 | 8.98 / 10.3 / 0.46 | 9.2 / 85.9 / 7.3 | 134 / 240 / 28.7 | 6.86 / 10.8 / 0.97
150 | 21.7 / 23.7 / 0.7 | 2.36 / 4.86 / 1.22 | 9.08 / 10.4 / 0.44 | 10.8 / 86.6 / 8.4 | 146 / 308 / 34.2 | 8.41 / 17.1 / 1.46

###  IV-C Environment with obstacles

We add obstacles to the environment as well as delay to see how our planner
performs as the communication latency increases. The obstacles have already
been inflated by the agent’s radius at their generation. We test for $8$ and
$12$ agents. Furthermore, we change the diameter of each agent to $0.3$ m,
$v_{\text{samp}}=3.5$ m/s, $a_{\text{max}}=30$ m/s2, $j_{\text{max}}=60$ m/s3,
$N=7$ and $d_{\text{thresh}}=0.2$ m for experimental diversity. We generate
$70$ obstacles of size $0.2\times 0.2\times 1.5$ m with random positions at
each simulation run (uniform distribution - Fig. 6(b), 6(c)). We do 10
simulation runs for each latency $dt=0,50,100$, and $150$ ms. The performance
metrics used are the distance traversed by each agent, the flight velocity and
time, the computation time, and the acceleration and jerk costs. The mean /
max / standard deviation of each metric are shown in Tab. III.

In all test runs for 8 and 12 agents, all agents were able to reach their
intended goal/destination safely i.e. the safety distance between the agents
was not violated and they did not get stuck along the way.

For 8 agents, the results for $dt=0$ ms and $dt=50$ ms are similar. This is
due to the fact that in both cases, all agents receive the trajectories before
the start of the next planning iteration since the maximum computation time is
below $50$ ms. The results for $dt=100$ ms and $dt=150$ ms are also similar
due to the same reason: in both cases, all agents receive the trajectories of
other agents every 2 planning iterations (the planning period is effectively
$2h$ due to our latency handling algorithm 1).

For 8 and 12 agents, the jerk cost and computation time both increase as the
latency increases. This is due to the more frequent slowdown of each agent as
the latency increases. The slowdown is due to passing through narrow spaces
and avoiding other agents at the same time as well as the latency handling
mechanism (see video link after the abstract).

##  V Conclusions and Future Works

In this paper, we presented an improved decentralized, real-time, and
synchronous framework for multi-agent planning. The method improves on our
previous work [1] by making it fully online and suitable for real-world
applications (the global path planning and Safe Corridor generation steps were
done offline in [1]). Furthermore, we added a mechanism to handle arbitrary
communication latency and adapt the planning frequency accordingly. Our
previous work was only able to handle communication latency when it is lower
than a predetermined threshold. We compared our work to 3 state-of-the-art
multi-agent planning methods: Ego-Swarm [8], MADER [7] and RMADER [13]. We
showed that our planner generates the safest trajectories with a $0$%
collision rate. Furthermore, it generates smoother and faster trajectories
than the only other safe and latency robust planner (RMADER) while also being
at least $3\times$ more computationally efficient.

In the future, we plan on implementing our planning method on embedded drone
systems for swarm autonomous navigation. This would require implementing
relative localization algorithms between agents, obstacle detection for
collision avoidance, as well as a communication mechanism for broadcasting
information between agents. Finally, we intend on developing a formation
flight version of our planner. This can be done by adding a cost to the
objective function of our planner that makes agents preserve a predefined
shape.

## References

  * [1] C. Toumieh and A. Lambert, “Decentralized multi-agent planning using model predictive control and time-aware safe corridors,” _IEEE Robotics and Automation Letters_ , pp. 1–8, 2022.
  * [2] W. Hönig, J. A. Preiss, T. K. S. Kumar, G. S. Sukhatme, and N. Ayanian, “Trajectory planning for quadrotor swarms,” _IEEE Transactions on Robotics_ , vol. 34, no. 4, pp. 856–869, 2018.
  * [3] J. Park, D. Kim, G. C. Kim, D. Oh, and H. J. Kim, “Online distributed trajectory planning for quadrotor swarm with feasibility guarantee using linear safe corridor,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 4869–4876, 2022.
  * [4] H. Zhu and J. Alonso-Mora, “B-uavc: Buffered uncertainty-aware voronoi cells for probabilistic multi-robot collision avoidance,” in _2019 International Symposium on Multi-Robot and Multi-Agent Systems (MRS)_ , 2019, pp. 162–168.
  * [5] D. Zhou, Z. Wang, S. Bandyopadhyay, and M. Schwager, “Fast, on-line collision avoidance for dynamic vehicles using buffered voronoi cells,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 2, pp. 1047–1054, 2017.
  * [6] C. E. Luis, M. Vukosavljev, and A. P. Schoellig, “Online trajectory generation with distributed model predictive control for multi-robot motion planning,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 2, pp. 604–611, 2020.
  * [7] J. Tordesillas and J. P. How, “Mader: Trajectory planner in multiagent and dynamic environments,” _IEEE Transactions on Robotics_ , pp. 1–14, 2021\.
  * [8] X. Zhou, J. Zhu, H. Zhou, C. Xu, and F. Gao, “Ego-swarm: A fully autonomous and decentralized quadrotor swarm system in cluttered environments,” _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 4101–4107, 2021.
  * [9] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3d mapping framework based on octrees,” _Autonomous robots_ , vol. 34, no. 3, pp. 189–206, 2013.
  * [10] E. Soria, F. Schiano, and D. Floreano, “Distributed predictive drone swarms in cluttered environments,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 1, pp. 73–80, 2021.
  * [11] C. Toumieh and A. Lambert, “Gpu accelerated voxel grid generation for fast mav exploration,” _under review for The Journal of Intelligent and Robotic Systems_ , 2021.
  * [12] B. Senbaslar and G. Sukhatme, “Asynchronous real-time decentralized multi-robot trajectory planning,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2022)_ , 2022.
  * [13] K. Kondo, J. Tordesillas, R. Figueroa, J. Rached, J. Merkel, P. C. Lusk, and J. P. How, “Robust mader: Decentralized and asynchronous multiagent trajectory planner robust to communication delay,” _arXiv preprint arXiv:2209.13667_ , 2022.
  * [14] C. Toumieh and A. Lambert, “High-speed planning in unknown environments for multirotors considering drag,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 7844–7850.
  * [15] D. D. Harabor and A. Grastien, “Online graph pruning for pathfinding on grid maps,” in _Twenty-Fifth AAAI Conference on Artificial Intelligence_ , 2011\.
  * [16] C. Toumieh and A. Lambert, “Near time-optimal trajectory generation for multirotors using numerical optimization and safe corridors,” _Journal of Intelligent & Robotic Systems_, vol. 105, no. 1, pp. 1–10, 2022.
  * [17] R. Deits and R. Tedrake, “Computing large convex regions of obstacle-free space through semidefinite programming,” in _Algorithmic foundations of robotics XI_. Springer, 2015, pp. 109–124.
  * [18] S. Liu, M. Watterson, K. Mohta, K. Sun, S. Bhattacharya, C. J. Taylor, and V. Kumar, “Planning dynamically feasible trajectories for quadrotors using safe flight corridors in 3-d complex environments,” _IEEE Robotics and Automation Letters_ , vol. 2, no. 3, pp. 1688–1695, 2017.
  * [19] C. Toumieh and A. Lambert, “Voxel-grid based convex decomposition of 3d space for safe corridor generation,” _Journal of Intelligent & Robotic Systems_, vol. 105, no. 4, pp. 1–13, 2022.
  * [20] ——, “Shape-aware safe corridors generation using voxel grids,” _arXiv e-prints_ , pp. arXiv–2208, 2022.
  * [21] J. Park, I. Jang, and H. J. Kim, “Decentralized deadlock-free trajectory planning for quadrotor swarm in obstacle-rich environments - extended version,” _ArXiv_ , vol. abs/2209.09447, 2022.