Counterfactual fairness for small subgroups–References

# Counterfactual fairness for small subgroups

Solvejg Wastvedt∗<EMAIL_ADDRESS>Jared D Huling  and Julian Wolfson
Division of Biostatistics  University of Minnesota  Minneapolis  Minnesota
U.S.A

###### Abstract

While methods for measuring and correcting differential performance in risk
prediction models have proliferated in recent years, most existing techniques
can only be used to assess fairness across relatively large subgroups. The
purpose of algorithmic fairness efforts is often to redress discrimination
against groups that are both marginalized and small, so this sample size
limitation often prevents existing techniques from accomplishing their main
aim. We take a three-pronged approach to address the problem of quantifying
fairness with small subgroups. First, we propose new estimands built on the
“counterfactual fairness” framework that leverage information across groups.
Second, we estimate these quantities using a larger volume of data than
existing techniques. Finally, we propose a novel data borrowing approach to
incorporate “external data” that lacks outcomes and predictions but contains
covariate and group membership information. This less stringent requirement on
the external data allows for more possibilities for external data sources. We
demonstrate practical application of our estimators to a risk prediction model
used by a major Midwestern health system during the COVID-19 pandemic.

###### keywords:

Algorithmic fairness, causal inference, risk prediction, small subgroups.

##  1 Introduction

Increasing use of complex risk prediction models in health care settings has
drawn attention to both the opportunities and the challenges such models
present. Clinical risk prediction models can improve care for patients and
efficiency for providers by personalizing treatment, identifying high-risk
patients for early intervention, and more. However, the models, which often
use opaque techniques that are not readily understandable by patients and
providers, have also demonstrated the potential to create and entrench health
inequities. One well-documented strand of this problem is the potential for
risk prediction models to perform more poorly for some groups than others. In
clinical settings, inaccurate model predictions can lead to sub-optimal
assignment of treatment and worse health outcomes for the affected groups,
which are often those already marginalized in society (Obermeyer et al.,
2019).

While techniques for assessing and correcting model bias have proliferated in
recent years (e.g., Castelnovo et al. (2022), Chen et al. (2023)), much of
this work does not address a major challenge in the clinical setting: limited
sample size in the smallest groups where model performance is to be assessed
to ascertain fairness. As an example, our previous work, which proposes
estimators of intersectional unfairness adapted for clinical situations in
which a treatment is in use, requires a reasonably large sample size in all
groups to obtain useful precision (Wastvedt et al., 2023). In that work, we
analyzed a risk prediction model used during the COVID-19 pandemic to help
determine whether a patient should be transferred to a COVID-19 cohort
hospital. We looked at unfairness across intersecting race and age categories,
but we were only able to separate race into white, Black or African American,
and a group for all other races because of limited sample size. While this
analysis reveals differing model error rates by age and for patients
identifying as white vs. Black or African American, the “all other races
groups” has little utility. The experiences of individuals in this group,
which includes for example both Asian and American Indian patients, differ
widely. It is likely that the model performs differently for the subgroups
encompassed in this category, but those error rate differences disappear when
the subgroups are examined together. More broadly, lumping small subgroups
together in this fashion can erase the experiences of these patients from
analyses. The main purpose of fairness assessment is to uncover mistreatment
of marginalized groups which are often, though not always, numeric minorities
as well. Thus it is crucial to the goals of fairness methods that they
prioritize the ability to obtain results for small subgroups.

In this work, we address the problem of limited sample size by proposing new
estimators that are estimable and have reasonable precision even for very
small groups. We extend the COVID-19 model application presented in our
earlier work to show the utility of these new estimators in situations where
existing methods fail. The COVID-19 risk prediction model we analyze is one of
many that were used to help guide allocation of resources in high-pressure
scenarios as case numbers rose. The model we analyze was trained by the health
system on 1,469 adults who tested positive for SARS-CoV-2. Features used for
prediction included patient demographic variables, home medications, and
medical conditions. Data available for our after-the-fact assessment of this
model is limited because of the ever-changing nature of COVID-19 treatment and
the relatively short period, during the height of the pandemic, over which the
model was used. As is common in clinical risk prediction scenarios, the
population for which the model was used is small in comparison to the
population served by the health system. Thus there is potential to borrow
information from other available data drawn from the same population in which
risk scores and outcomes are not recorded; we explore this possibility in our
proposed methods.

In the sections that follow we first describe the _counterfactual error rates_
, model performance metrics developed in existing work (Coston, Mishler,
Kennedy, and Chouldechova, 2020; Mishler, Kennedy, and Chouldechova, 2021)
that are well-suited to clinical applications. We then take a three-pronged
approach to improving these metrics’ performance for small subgroups. First,
we re-formulate the error rates to leverage information across groups by
incorporating the overall error rate. We then use the causal identification
process to arrive at estimators that use more of the data than existing
methods. Finally, we propose a novel data borrowing procedure that uses
“external data”, such as that available in our COVID-19 risk model
application, to further improve estimation.

##  2 Statistical framework

Following algorithmic fairness convention, define a _protected characteristic_
as any grouping, such as race or gender, along which we wish to measure
discrimination. Let $\boldsymbol{A}$ denote a vector of categorical protected
characteristics, in which each element indicates group membership for one
characteristic. Let $m$ be the number of protected characteristics that we
wish to consider. Denote the characteristics $A_{j}$, $j\in\\{1,...,m\\}$ and
assume each $A_{j}$ is a categorical variable with a finite number of levels,
the set of which is denoted $\mathcal{A}_{j}$. Let
$\boldsymbol{A}=\\{A_{1},A_{2},...,A_{m}\\}^{T}\in\mathcal{A}$ contain all
protected characteristics of interest, where $\mathcal{A}$ is the set of all
possible combinations of all levels of the $m$ characteristics.

To finish notation, let $S$ denote a binary risk prediction. Although many
clinical risk models produce predicted probabilities, we are using the binary
prediction derived by selecting a cut-off threshold. Let $D$ denote a binary
treatment assignment and $Y$ a binary outcome such as an adverse health event.
Under a binary treatment, there are two potential outcomes: $Y^{0}$, the
outcome under no treatment, and $Y^{1}$, the outcome under treatment. We focus
on the $Y^{0}$ outcome as it represents patients’ baseline risk and thereby is
informative for guiding treatment.

Existing work (Mishler et al., 2021) defines the following counterfactual
versions of the common false positive rate and false negative rate model
performance metrics for a single protected characteristic. In our previous
work, we apply the definitions to the vector of protected characteristics
$\boldsymbol{A}=\boldsymbol{a}$.

###### Definition 2.1

The _counterfactual false positive rate_ of a predictor $S$ for the group
having protected characteristic vector $\boldsymbol{A}=\boldsymbol{a}$,
denoted $cFPR(S,\boldsymbol{a})$, is equal to
$Pr(S=1|Y^{0}=0,\boldsymbol{A}=\boldsymbol{a})$. The _counterfactual false
negative rate_ , $cFNR(S,\boldsymbol{a})$, is equal to
$Pr(S=0|Y^{0}=1,\boldsymbol{A}=\boldsymbol{a})$.

Our previous work proposed the following weighted estimators of
$cFPR(S,\boldsymbol{a})$ and $cFNR(S,\boldsymbol{a})$:

$\displaystyle\widehat{cFPR}(S,\boldsymbol{a})$
$\displaystyle=\frac{\sum_{i=1}^{n}[(1-D_{i})S_{i}(1-Y_{i})\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})/(1-\hat{\pi}_{i})]}{\sum_{i=1}^{n}[(1-D_{i})(1-Y_{i})\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})/(1-\hat{\pi}_{i})]}$
(1) $\displaystyle\widehat{cFNR}(S,\boldsymbol{a})$
$\displaystyle=\frac{\sum_{i=1}^{n}[(1-D_{i})(1-S_{i})Y_{i}\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})/(1-\hat{\pi}_{i})]}{\sum_{i=1}^{n}[(1-D_{i})Y_{i}\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})/(1-\hat{\pi}_{i})]}$
(2)

.

The _counterfactual error rate differences_ , denoted $\Delta^{+}(S)$ and
$\Delta^{-}(S)$, quantify the differences in counterfactual false positive and
false negative rates between two protected characteristic groups.

##  3 New estimands to borrow strength across groups

The estimators for $cFPR(S,\boldsymbol{a})$ and $cFNR(S,\boldsymbol{a})$ in
equations (1) and (2) present a challenge in practical applications with
limited sample size. The estimators are weighted averages over a subset of a
subgroup defined by intersecting protected characteristics. In the case of
$cFPR(S,\boldsymbol{a})$, only observations with $D_{i}=0$, $S_{i}=1$, and
$Y_{i}=0$ contribute non-zero components to the average. The estimator for
$cFNR(S,\boldsymbol{a})$ is analogous. Restricting the estimation to this
slice of what may already be a small protected subgroup can make the
estimators unstable.

To address this issue, we propose alternative estimators for the
counterfactual error rates that incorporate more of the data. Our three-
pronged approach begins by reformulating the estimands for group-specific
error rates to leverage the overall error rate. Next, we make an additional
assumption during causal identification to reduce reliance on indicator
variables. Finally, the next section explains potential use of an auxiliary
data set during estimation to further improve performance.

Our re-formulation of the counterfactual error rates draws on Efron (2010),
who propose rewriting the false discovery rate (FDR) of a predictor $S$ for a
particular subgroup as a ratio of probabilities multiplied by the overall FDR.
Following a similar process, we obtain alternate expressions for the
counterfactual error rates given in Proposition 3.1 (proof in Web Appendix A).

###### Proposition 3.1

Given two binary variables $S$ and $Y^{0}$ and a subgroup vector
$\boldsymbol{A}=\boldsymbol{a}$, the counterfactual false negative rate and
counterfactual false positive rate of $S$ in group $\boldsymbol{a}$ can be
rewritten as the following:

$\displaystyle cFPR(S,\boldsymbol{a})$
$\displaystyle=cFPR(S)\frac{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=0,S=1)}{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=0)}$
(3) $\displaystyle cFNR(S,\boldsymbol{a})$
$\displaystyle=cFNR(S)\frac{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=1,S=0)}{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=1)}$
(4)

.

We then identify our estimands from observed data using standard causal
inference assumptions plus one additional assumption that allows us to avoid
many of the indicator variables used in Wastvedt et al. (2023), thereby
employing a larger volume of data. We assume $Y^{0}$ is independent of
$\boldsymbol{A}$ given $X$, or that we have collected sufficient covariates
such that protected group membership gives no additional information about the
probability of $Y^{0}$. A full list of assumptions and simulations
demonstrating that our methods are robust to violations of the additional
assumption are in Web Appendices B and D.

Using these assumptions, we break down the ratio of conditional probabilities
in each error rate expression in Proposition 3.2 (proof in Web Appendix B).
Let $X$ be a vector of observed covariates and define the propensity score
function for a given protected group as $\pi=P(D=1|A,X,S)$.

###### Proposition 3.2

The conditional probability ratios in expressions (3) and (4) are identified
as follows for $Y^{0}=y\in\\{0,1\\}$ and $S=s\in\\{0,1\\}$:

$\frac{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=y,S=s)}{P(\boldsymbol{A}=\boldsymbol{a}|Y^{0}=y)}=\frac{E[\mu_{0}(y,s,X)\mathbbm{1}(\boldsymbol{A}=\boldsymbol{a})\mathbbm{1}(S=s)]/E[\mu_{0}(y,s,X)\mathbbm{1}(S=s)]}{E[\mu_{0}^{*}(y,X)h(X,\boldsymbol{a})]/E[\mu_{0}^{*}(y,X)]}\text{,}$

where $\mu_{0}(y,s,X)=P(Y=y|D=0,S=s,X)$ and $\mu_{0}^{*}(y,X)=P(Y=y|D=0,X)$.
The $h$ function is defined as
$h(x,\boldsymbol{a})=P(\boldsymbol{A}=\boldsymbol{a}|X=x)$.

While the estimators in Wastvedt et al. (2023) limited data to specific
intersections of $D$, $S$, $Y$, and $\boldsymbol{A}$, Proposition 3.1 allows
us to estimate the denominator using all the data and only restrict by the
protected characteristic vector ($\boldsymbol{A}$) and model prediction ($S$)
in the numerator. Drawing from more data reduces our estimators’ variance and
enables estimation even in very small subgroups.

The overall counterfactual error rate components of equations (3) and (4) use
the identification result in Wastvedt et al. (2023), generalized to include
all observations:
$cFPR(S)=\frac{E[(1-D)S(1-Y)/(1-\pi(\boldsymbol{A},X,S))]}{E[(1-D)(1-Y)/(1-\pi(\boldsymbol{A},X,S))]}$;
$cFNR(S)=\frac{E[(1-D)(1-S)Y/(1-\pi(\boldsymbol{A},X,S))]}{E[(1-D)Y/(1-\pi(\boldsymbol{A},X,S))]}$.

##  4 Estimation and data borrowing

In this section, we demonstrate use of an auxiliary, external data set to aid
in estimation of $h(X,\boldsymbol{a})=P(\boldsymbol{A}=\boldsymbol{a}|X)$.
This procedure can reduce variance and bias compared to estimation using only
the main data. We then propose regression estimators for the remaining
components of the probability ratios in Proposition 3.2.

Throughout, let $\\{\boldsymbol{A}_{i},D_{i},Y_{i},X_{i},S_{i}\\}$,
$i=1,...,n$ be the observed data and binary risk predictions. Then the
counterfactual error rates for protected group $\boldsymbol{A}=\boldsymbol{a}$
are estimated as follows, where for conciseness we assume that $\hat{\mu}_{0}$
and $\hat{\mu}_{0}^{*}$ are the estimated versions with $y=1$:

$\displaystyle\widehat{cFPR(S,\boldsymbol{a})}$
$\displaystyle=\widehat{cFPR(S)}\frac{\sum_{i=1}^{n}(1-\hat{\mu}_{0}(s=1,X_{i}))\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})S_{i}/\sum_{i=1}^{n}(1-\hat{\mu}_{0}(s=1,X_{i}))S_{i}}{\sum_{i=1}^{n}(1-\hat{\mu}_{0}^{*}(X_{i}))\hat{h}(X_{i},\boldsymbol{a})/\sum_{i=1}^{n}(1-\hat{\mu}_{0}^{*}(X_{i}))}$
(5) $\displaystyle\widehat{cFNR(S,\boldsymbol{a})}$
$\displaystyle=\widehat{cFNR(S)}\frac{\sum_{i=1}^{n}\hat{\mu}_{0}(s=0,X_{i})\mathbbm{1}(\boldsymbol{A}_{i}=\boldsymbol{a})(1-S_{i})/\sum_{i=1}^{n}\hat{\mu}_{0}(s=0,X_{i})(1-S_{i})}{\sum_{i=1}^{n}\hat{\mu}_{0}^{*}(X_{i})\hat{h}(X_{i},\boldsymbol{a})/\sum_{i=1}^{n}\hat{\mu}_{0}^{*}(X_{i})}$
(6)

.

###  4.1 Data borrowing to estimate group membership probabilities

Estimating the group membership probabilities $\hat{h}(X,\boldsymbol{a})$
involves neither prediction ($S$) nor outcome ($Y$) information. In many
clinical risk prediction situations, large volumes of data are available in
which the outcome and/or prediction are not present, but a rich set of
covariates are. For example, in a health system, this external data could be
patient records in the same electronic health record system who were not
screened for the adverse event or who did not receive risk predictions.
Because of the larger sample size, this external data will typically have
better representation from small subgroups. If the distribution of
$P(\boldsymbol{A}=\boldsymbol{a}|X)$ in the external data is similar enough to
that in the test (or “internal”) data, we can use the external data to help
estimate group membership probabilities.

To quantify “similar enough”, we employ an adaptive data borrowing method that
borrows more information if the distributions are more closely aligned and
less (or none) if they differ. Let $\hat{h}_{E}$ be the internal data group
membership probabilities estimated using a model trained on the external data.
Let $\hat{h}_{I}$ be the probabilities fit and predicted using only the
internal data. We maximize the predictive performance of the function
$\hat{h}^{*}=\hat{\alpha}\hat{h}_{E}+(1-\hat{\alpha})\hat{h}_{I}$ to find the
borrowing parameter $\hat{\alpha}$. Various metrics could be used for
predictive performance. We use the Brier score because it encompasses both
concordance and calibration; however, we achieve similar results in
simulations using multi-class area under the ROC curve (AUC). The estimates
$\hat{h}^{*}_{i}$ are then used in equations (5) and (6).

###  4.2 Estimation of remaining quantities

The estimators in equations (5) and (6) require several additional parameters.
The two outcome models, $\mu_{0}(X_{i},S=1)$ and $\mu^{*}_{0}(X_{i})$, can be
estimated using logistic regression or a more flexible method suitable for
binary data. The overall counterfactual error rates are estimated using the
weighted estimators proposed in Wastvedt et al. (2023):
$\widehat{cFPR(S)}=\frac{\sum_{i=1}^{n}[(1-D_{i})S_{i}(1-Y_{i})/(1-\hat{\pi}_{i})]}{\sum_{i=1}^{n}[(1-D_{i})(1-Y_{i})/(1-\hat{\pi}_{i})]}$,
$\widehat{cFNR(S)}=\frac{\sum_{i=1}^{n}[(1-D_{i})(1-S_{i})Y_{i}/(1-\hat{\pi}_{i})]}{\sum_{i=1}^{n}[(1-D_{i})Y_{i}/(1-\hat{\pi}_{i})]}$.
The propensity score, $\pi(\boldsymbol{A},X,S)$, can be estimated using
regression methods similar to the outcome models.

Below we compare estimation of propensity score and outcome models using a
generalized linear model and a Super Learner (Laan et al., 2007), a more
complex ensemble approach that combines multiple machine learning models.
Because of the simplicity of the generalized linear models, we use the full
sample without any sample splitting or cross-fitting. For the ensemble models,
we use a 10-fold cross-fitting approach. For each fold, we fit the three
nuisance parameter models on the remaining data to obtain predictions for the
held-out fold.

##  5 Simulations

In this section we demonstrate the benefits of our proposed estimators
compared to existing methods. We compare to the counterfactual estimators in
(Wastvedt et al., 2023) that are not adapted for small subgroups.

All simulations consider two binary protected characteristics, $A_{1}$ and
$A_{2}$, where for each level $1$ is less common. Thus the group $A_{1}=0$,
$A_{2}=0$ is the numeric majority, and the group $A_{1}=1$, $A_{2}=1$ is the
numeric minority. We denote the other groups $M1$ and $M2$ by the protected
characteristic ($A_{1}$ or $A_{2}$) that is equal to $1$. Overall internal
data sample size is limited in all scenarios to create small groups
($N_{internal}=\\{50,100,150,200\\})$. External data is large
($N_{external}=10,000$) such that all groups are adequately represented in the
external sample. We focus on $cFNR(S,\boldsymbol{a})$ for conciseness; our
conclusions also hold for $cFPR(S,\boldsymbol{a})$.

Following the simulation framework in Mishler et al. (2021), we first generate
a set of data ($N=1,000)$ for training a random forest risk prediction model.
We generate a validation data set ($N=50,000$) to determine the true error
rates of the risk model. We use the same risk prediction model for all
scenarios. For each scenario, we generate internal and external data sets and
estimate our unfairness metrics, considering multiclass area under the ROC
curve (AUC) and Brier score as performance metrics for data borrowing. In all
scenarios, the internal and external $h(X,\boldsymbol{a})$ models are neural
networks with a single hidden layer, $100$ units, and a weight decay parameter
of $1$ (R package nnet). Full data generation details are in Web Appendix C.

###  5.1 Variance reduction for small subgroups

This section compares our proposed and existing estimators under internal
sample sizes that are at the lower end of what would typically be encountered
in a clinical risk prediction assessment setting
($N_{int}\in\\{100,200,500,1000,2000\\}$). We use the Brier score for the data
borrowing metric. Outcome models (proposed estimators) and the propensity
score model (comparison estimators) are correctly specified GLMs.

Figure 1 shows notable variance reduction with our new estimators across
sample sizes. Error bars without a mean (shape) indicate the presence of NAs
in the replications, i.e. insufficient data for estimation. This occurs with
our comparison estimators in all groups at $N_{int}=200$ and in non-majority
groups at higher $N_{int}$. In contrast, both our proposed internal and
borrowing estimators are able to be calculated for all groups at all sample
sizes. Comparing internal and borrowing shows little gain with data borrowing
under these particular scenarios.

###  5.2 Data borrowing adapts to external data agreement

This section investigates our data borrowing procedure’s ability to respond to
varying levels of disagreement in internal and external data distributions. In
this section we fix $N_{int}=1000$ and manipulate the level of agreement
between internal and external data by multiplying the external data
coefficients for $P(\boldsymbol{A}=\boldsymbol{a}|X)$ by a constant in
$[-1,1]$.

Figure 2(a) shows that $\alpha$ decreases as $b$ moves away from $1$. Recall
that a lower $\alpha$ means less weight on the external data predictions,
i.e., less borrowing. Figure 2(b) focuses on the minority group $cFNR$ to show
that across all distributions, our data borrowing estimators have bias that is
less than or equal to the internal data estimation estimator and reasonable
variance.

###  5.3 Benefit of data borrowing under complex scenarios

Next we examine scenarios where data borrowing provides gains in bias and
variance reduction compared to the internal data version of our estimators.
Such gains occur when the internal $h(X,\boldsymbol{a})$ model is difficult to
estimate. We use a more complex data generation scenario with an increasing
number of noise variables added to the $10$ informative components of $X$
before generating group membership, $P(\boldsymbol{A}=\boldsymbol{a}|X)$, the
$Y^{0}$ potential outcome, $P(Y^{0}=1|X,\boldsymbol{A})$, and the treatment,
$P(D=1|X,\boldsymbol{A},S)$. We consider these scenarios with and without 2-
and 3-way interactions among four of the $X$ used to generate these
probabilities. For simplicity, we focus on the minority group and Brier score
for the data borrowing metric. The internal sample size is $500$ throughout.

Figure 3 shows that both with and without $X$ interactions, bias is low and
consistent across increasing $X$ noise up to 20 $X$ components. Data borrowing
does not notably affect bias except a small reduction under the no-
interaction, 40 $X$ components scenario.

##  6 COVID-19 risk prediction fairness

During the height of the COVID-19 pandemic, risk prediction models were
commonly used to help health systems allocate scarce resources. In addition to
their potential benefits, these models raised the challenge of ensuring model
performance, and thus resource allocation, was equitable. We studied one such
model at a major Midwestern health system that was used to help determine
whether a patient should be transferred to a COVID-19 cohort hospital. In our
previous work, we analyzed this model’s counterfactual error rates across the
intersection of age and race, grouping patient self-reported race into white,
Black or African American, and all other races. As we noted in that paper,
grouping all patients who identify as neither white nor Black/African American
together severely limits the utility of the analysis. Here we show how our
proposed estimators allow us to disaggregate race and use an additional
ethnicity variable to gain a more nuanced picture of model performance. We
focus on a single protected characteristic to demonstrate the wide
applicability of our proposed methods; our estimators can also be used with
multiple, intersecting characteristics.

The risk prediction model was trained on 1,469 adult patients with confirmed
or symptomatic suspected COVID-19. Our “internal” data for evaluating the
model consisted of 3,649 adult SARS-CoV-2 patients from the same health system
who tested positive between 10/27/2020 and 1/9/2022. Deploying our new
estimators allowed us to group patient self-reported race and ethnicity
variables into five categories: Hispanic or Latino (6.3%); and non-
Hispanic/Latino American Indian or Alaska Native (1.5%), Asian (4.2%), Black
or African American (13.9%), and white (74.1%). We removed patients who marked
two or more races, “other” or “all other races” since these records lacked
sufficient information to group race/ethnicity. We chose this schema to
balance disaggregation with a large enough sample size to obtain useful
estimates. Our race categories align with U.S. Census Bureau categories with
one exception: we did not have sufficient sample size to break out the Native
Hawaiian or Other Pacific Islander category since it had 14 total
observations, so we grouped these patients with the Asian category. Likewise,
cross-cutting the data by both race and Hispanic or Latino ethnicity was not
possible given our sample size, so we opted to consider Hispanic or Latino
ethnicity by disaggregating it from the racial groups. We note that the data
used in this application does offer relatively rich detail on patients’ self-
reported race and ethnicity, with 8 and 48 unique responses to the race and
ethnicity questions, respectively. Future applications using data sources like
this with larger samples could take fuller advantage of the available detail
by considering race and ethnicity separately and with more categories.

We chose a cutoff of $0.15$, approximately the $80^{th}$ percentile, for
dichotomizing the risk score and used 30-day inpatient readmission or
mortality as our outcome and transfer to the cohort hospital as our treatment
variable. Covariates comprised comorbidities, home medications, number of
prior emergency department visits, and labs and vitals. We excluded covariates
with greater than 2/3 missing and performed random forest imputation
(missForest package). We considered the same set of covariates for propensity
score, outcome, and $h(X,\boldsymbol{a})$ modeling, performing lasso selection
separately for each model to select the most relevant variables. As in our
prior work, we used a logistic regression propensity score model. Our outcome
models were also logistic regressions, and the internal $h(X,\boldsymbol{a})$
model was a single-layer neural network with 50 units (package nnet).

Our external data comprised 8,449 patients from the same health system without
a matching SARS-CoV-2 positive test and risk model score. External data
patients were grouped into the same five race/ethnicity categories: Hispanic
or Latino (8.6%); and non-Hispanic/Latino American Indian or Alaska Native
(1%), Asian (11.2%), Black or African American (20%), and white (59.2%). While
not all internal data covariates were available in the external data,
covariates comprised the same general categories. We excluded covariates with
greater than 2/3 missing and performed random forest imputation. We performed
lasso selection to choose the most important covariates for the external
$h(X,\boldsymbol{a})$ model and then fit a single-layer neural network with
100 units using the selected covariates. We used the Brier score data
borrowing metric which gave $\alpha=0.004$. We obtained standard errors using
the rescaled bootstrap method proposed in Wastvedt et al. (2023) and
confidence intervals using t-intervals truncated at $0$ and $1$. All analyses
were done in R (version 4.2.3, R Core Team 2023).

Our proposed estimators enable estimation for small subgroups where estimation
failed using the comparison estimators. Figure 4 shows that internal and
borrowing estimators reduced variance substantially for the Black or African
American group and slightly for the white group. The smallest group, American
Indian or Alaska Native, was too small to obtain a confidence interval and
obtained the unlikely point estimate of $1$ using the comparison method. Our
new estimators obtain an interval and a more reasonable point estimate,
although the interval covers the entire $(0,1)$ range and thus has minimal
utility. This wide interval demonstrates the sample size limitations that
still occur using our proposed estimators.

The Asian and Hispanic or Latino groups had point estimates that differ
substantially under the new and comparison methods, which may be due to
undercoverage issues for the comparison method. Estimates for these groups
using the comparison method are also highly variable because of extremely
small counts ($<10$) in the false negative cell of the confusion matrix, which
is the relevant cell for the comparison $cFNR$ estimate. Note also that while
the confidence interval for the Asian group appears smaller under the
comparison method, this is due to the truncation of intervals at zero.
Intervals affected by truncation are noted in the caption to Figure 4.

##  7 Discussion

In this paper we took a three-pronged approach to addressing the common
challenge of small subgroups in risk prediction model fairness assessments.
First, we proposed new estimands that reformulate the counterfactual error
rates to borrow strength across groups. We showed how these estimands can be
identified from observed data in a manner that draws on more of the data than
existing alternatives. Finally, we proposed a novel data borrowing method with
the potential to leverage external data to further improve estimation.

Our methods have limitations: most importantly, subgroups can still be too
small for the methods to produce useful estimates. When sample size is too
small, confidence intervals may be too wide, as in the American Indian or
Alaska Native group in our application (Section 6), or estimates may differ
substantially between our proposed method and the comparison method. However,
this sample size cutoff is fairly small, as demonstrated by the useful
estimates produced for the Black or African American group in our application,
a group which only had $11$ observations in the false negative cell of the
confusion matrix (see Web Table 1 for confusion matrices by group). We refrain
from offering a specific sample size cut-off for our methods since this number
will depend on the richness of the covariates available in a given
application, both for model fitting and for meeting our additional assumption
described in Section 3. Our proposed estimators also require practitioners to
fit more models than existing alternatives which is both more time-consuming
and provides more opportunities for model mis-specification. Nevertheless, we
find little impact of model mis-specification under the scenarios considered
in our simulations (Web Appendix D).

Finally, our estimators allow practitioners to look at more disaggregated
protected groups, which we hope will facilitate more nuanced consideration in
fairness assessments of groups that are both marginalized and small. In the
case of race and ethnicity, as in our COVID-19 application, our estimators are
an improvement over those that require grouping disparate subgroups or
excluding small groups altogether. However, given the prevalence of using
racial categories as a protected characteristic in this field, it is important
to recognize the limitations inherent in any race-based fairness assessment.
No matter how disaggregated, racial categories are not intrinsic identities;
rather, they are socially constructed labels that comprise “a system of
inherently unequal status categories” (Benthall and Haynes, 2019). In a
critique of traditional approaches to algorithmic fairness, Weinberg (2022)
frames the problem as one of abstraction: by treating race as an intrinsic
quality, researchers minimize the structural factors that create and maintain
racial inequality. Whether conducted with our proposed metrics or others, a
fairness assessment involving race bears the responsibility to drive toward
undermining the systems of hierarchy and power imbalances that cause
algorithmic unfairness. Hanna et al. (2020) suggest considering descriptive
analyses of model performance across the many dimensions of race such as self-
identified, other-identified, phenotype, and others. Differing patterns of
unfairness can provide more insight into how to proceed to ameliorate
inequities. Finally, Hanna et al. (2020) emphasize the importance of looking
beyond the algorithm to promote justice at all parts of the process, including
the data models are trained on, the choice of where and on whom to deploy
them, and the way model predictions are incorporated into human practices.

##  8 Figures and Tables

Figure 1: Estimation of $cFNR(S,\boldsymbol{a})$ with varying internal data
sample size. Shapes show the mean of $500$ replications and error bars show
$95\%$-tile intervals. Dotted horizontal lines show the true value as obtained
using the validation set. Error bars without a mean (shape) indicate
insufficient data for estimation in at least one replication. In contrast to
the comparison estimators, our new estimators achieve reductions in bias and
variance and are usable at sample sizes where the comparison estimators cannot
be calculated.

(a) $\alpha$ values at varying levels of agreement between external and
internal data. As the data distributions diverge (left side), less information
is borrowed from the external data ($\alpha$ approaches zero).

(b) Bias in estimation of minority group $cFNR$ at varying levels of data
agreement. Right facet shows bias for internal data estimation (no borrowing).

Figure 2: Results of data borrowing with varying levels of agreement between
internal and external data distributions. The x-axis shows the factor by which
internal data coefficients for generating $P(\boldsymbol{A}=\boldsymbol{a}|X)$
were multiplied to obtain the external data coefficients. Shapes show the mean
of 500 replications of the estimation procedure, error bars show 95%-tile
intervals, and colors denote data borrowing metric. Figure 3: Bias in
estimation of minority group $cFNR$ under increasing noise in $X$ with
interactions among four of the $X$ (right) and no interactions (left). Shapes
show the mean of 500 replications of the estimation, and error bars show
95%-tile intervals. With no $X$ interactions, bias under data borrowing is
similar to internal-only estimation until 40 $X$, where borrowing performs
slightly better. With $X$ interactions, bias under data borrowing is similar
to internal-only estimation under all scenarios tested. For $\alpha$ values
showing how much information is borrowed in each scenario, see Web Appendix D.
Figure 4: Estimation of group $cFNR$ for 5 racial and ethnic groups using our
proposed estimators (internal and borrowing) and our comparison estimator.
Shapes show the mean of 200 bootstrap replications, and error bars show 95%
bootstrap t-intervals. Our proposed estimators reduce variance for the Black
or African American group and enable estimation for the American Indian or
Alaska Native group, although the confidence interval for this group covers
the entire $(0,1)$ range. Comparison estimators for Asian, Black or African
American, and Hispanic or Latino groups are truncated at $0$. New estimators
(internal and data borrowing) for the American Indian or Alaska Native group
are truncated at $1$.

## References

  * Benthall and Haynes (2019) Benthall, S. and Haynes, B. D. (2019).  Racial categories in machine learning.  In Proceedings of the Conference on Fairness, Accountability, and Transparency, pages 289–298, Atlanta, Georgia, USA. Association for Computing Machinery.
  * Castelnovo et al. (2022) Castelnovo, A., Crupi, R., Greco, G., Regoli, D., Penco, I. G., and Cosentini, A. C. (2022).  A clarification of the nuances in the fairness metrics landscape.  Scientific Reports 12, 4209\.
  * Chen et al. (2023) Chen, Z., Zhang, J. M., Sarro, F., and Harman, M. (2023).  A Comprehensive Empirical Study of Bias Mitigation Methods for Machine Learning Classifiers.  ACM transactions on software engineering and methodology 32, 1–30.
  * Coston et al. (2020) Coston, A., Mishler, A., Kennedy, E. H., and Chouldechova, A. (2020).  Counterfactual risk assessments, evaluation, and fairness.  In Proceedings of the Conference on Fairness, Accountability, and Transparency, pages 582–593, Barcelona, Spain. Association for Computing Machinery.
  * Efron (2010) Efron, B. (2010).  Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction.  Institute of Mathematical Statistics Monographs. Cambridge University Press, Cambridge.
  * Hanna et al. (2020) Hanna, A., Denton, E., Smart, A., and Smith-Loud, J. (2020).  Towards a Critical Race Methodology in Algorithmic Fairness.  In Proceedings of the Conference on Fairness, Accountability, and Transparency, pages 501–512, Barcelona, Spain. Association for Computing Machinery.
  * Laan et al. (2007) Laan, M. J. v. d., Polley, E. C., and Hubbard, A. E. (2007).  Super Learner.  Statistical Applications in Genetics and Molecular Biology 6,.
  * Mishler et al. (2021) Mishler, A., Kennedy, E. H., and Chouldechova, A. (2021).  Fairness in Risk Assessment Instruments: Post-Processing to Achieve Counterfactual Equalized Odds.  In Proceedings of the Conference on Fairness, Accountability, and Transparency, pages 386–400, Virtual Event, Canada. Association for Computing Machinery.
  * Obermeyer et al. (2019) Obermeyer, Z., Powers, B., Vogeli, C., and Mullainathan, S. (2019).  Dissecting racial bias in an algorithm used to manage the health of populations.  Science 366, 447–453.
  * Wastvedt et al. (2023) Wastvedt, S., Huling, J. D., and Wolfson, J. (2023).  An intersectional framework for counterfactual fairness in risk prediction.  Biostatistics page kxad021.
  * Weinberg (2022) Weinberg, L. (2022).  Rethinking Fairness: An Interdisciplinary Survey of Critiques of Hegemonic ML Fairness Approaches.  Journal of Artificial Intelligence Research 74, 75–109.


# Towards biomimicry of a bat-style perching maneuver on structures: the
manipulation of inertial dynamics

Alireza Ramezani Alireza Ramezai is with the Department of Electrical and
Computer Engineering, Northeastern University, Boston, MA 02115
<EMAIL_ADDRESS>
###### Abstract

The flight characteristics of bats remarkably have been overlooked in aerial
drone designs. Unlike other animals, bats leverage the manipulation of
inertial dynamics to exhibit aerial flip turns when they perch. Inspired by
this unique maneuver, this work develops and uses a tiny robot called Harpoon
to demonstrate that the preparation for upside-down landing is possible
through: 1) reorientation towards the landing surface through zero-angular-
momentum turns and 2) reaching to the surface through shooting a detachable
landing gear. The closed-loop manipulations of inertial dynamics takes place
based on a symplectic description of the dynamical system (body and
appendage), which is known to exhibit an excellent geometric conservation
properties.

##  I Introduction

The key hypothesis this work tries to inspect is: ”The feasibility of
performing flip turns, a well known attribute of a bat-like landing maneuver,
through the manipulation of inertial dynamics while excluding aerodynamic
forces.” Bats (and birds) possess no energy-hungry motors widely used in
flying robots for thrust vectoring yet they are more capable than any of these
systems when agility and, of course, energy efficiency of flight are
concerned.

Flying vertebrates apply the combination of inertial dynamics and aerodynamics
manipulations to showcase extremely agile maneuvers. Unlike rotary- and fixed-
wing systems wherein aerodynamic surface (i.e., ailerons, rudders, propellers,
etc.) come with the sole purpose of aerodynamic force adjustment, the wings
(also called appendages) in birds and bats possess more sophisticated roles.
It is known that birds perform zero-angular-momentum turns by making
differential adjustments (e.g., collapsing armwings) in the inertial forces
led by one wing versus the other. Or, bats apply a similar mechanism to
perform sharp banking turns [1].

Figure 1: Harpoon-V.2 developed at Northeastern University to test the idea of
inertial dynamics manipulation for fast aerial robot landing on ceilings and
steel frames in residential spaces.

Among these maneuvers, landing (or perching), which flying vertebrates do it
in one way or another for a variety of reasons (e.g., transition to walking,
resting on a perch, hanging from the ceiling of a cave, etc.), is an
interesting maneuver to take inspiration from for aerial robot designs.
Perching birds rotate their wings so that the aerodynamic drag is increased by
creating a high-pressure region inside of the wings and a low-pressure region
behind the wings. This brings the wings to a stalled condition at which point
the generated lift is equal to zero and the animal falls naturally while
employing the legs as a landing gear.

Bats do it in a radically different way. After the self-created stalled
condition, they manifest an acrobatic heels-above-head maneuver that involves
catapulting the lower body in a similar way that a free style swimmer flip
turns.

Perching insects and birds have been the source of inspiration and bio-mimicry
of them has led to interesting robot designs in recent years. Remarkably, the
bio-mimicry of bat-like landing is overlooked mainly because not only the
aerodynamics adjustments are involved but also unique design provisions are
required to allow for the manipulation of inertial dynamics.

A part from the ordeals associated with the design and control of a robot that
can land similar to a biological bat a number of unique applications can be
identified for these systems. For instance, in scenarios wherein maintaining a
high vantage point for extended times is demanded (surveillance and
reconnaissance) and limited power budget does not allow to hover for extended
time periods hanging from elevated structures such as steel frames in
buildings can allow these systems not only safely accomplish their missions
but also harvest energy within the time period the system is the most
vulnerable.

In addition to the exciting applications, a bat-style landing maneuver is
extremely rich in dynamics and control and its characteristics are overlooked.
Much of attention has been paid to simpler dynamics such as hovering and
straight flight. While mathematical models of insect-style, rotary- and fixed-
wing robots of varying size and complexity are relatively well developed,
models of airborne, fluidic-based vertebrates locomotion remain largely open
due to the complex body articulation involved in their flight.

The mainstream school of thought inspired by insect flight has conceptualized
wing as a mass-less, rigid structure, which is nearly planar and translates –
as a whole or in two-three rigid parts – through space. In this view, wings
possess no inertial effect, are fast that yield two-time-scale dynamics,
permit quasi-static external force descriptions, and tractable dynamical
system. Unfortunately, these paradigms fail to provide insight into airborne,
vertebrate locomotion and an ingredient of a more complete and biologically
realistic model is missing, that is, the manipulation of inertial dynamics.
The manipulation of inertial dynamics is an under-appreciated aspect in
existing paradigms.

The objective of this work is to manifest closed-loop aerial body
reorientation and preparation for landing through the manipulation of inertial
dynamics using a tiny ballistic robot whose characteristics are carefully
scaled to match that of a small bat. It is shown that despite a number of
prohibitive restrictions, and lack of multi-thruster designs typically found
in quadrotors, extremely fast body reorientations and preparations for landing
through inertial dynamics manipulations not only is possible but also is an
effective biologically meaningful solution.

A brief overview of the system will be presented followed by a hybrid-model
description of the landing maneuver. A closed-loop controller will be designed
and the simulation and experimental results will be reported at the end
followed by final remarks and conclusion.

##  II A hybrid-model description of upside-down perching

A little reflection reveals, however, that a perching bat dynamics like other
complex animal behavior, which emerges from complex interactions between
neural and sensory-motor systems, can be described with tractable mathematical
models. The legitimate question that is posed here is: ”Can a hybrid model
capture the dynamics?” What encourages me to pose the question directly stems
from inspecting high-speed images of the landing maneuver of bats and that in
a time envelope not exceeding about one fifth of a second before landing
completion (touch-down) the appendages are fully retracted strongly suggesting
that the aerodynamic forces are less likely to be as important as other
forces. Based on a similar observation, [2] conducted a simulation-based
research and concluded that the inertial dynamics manipulation likely are the
dominant players in a landing bat.

It is assumed for two modes: mode (1) with aerodynamics dominance (inertial
forces are not negligible but are less effective) and mode (2) within which
the inertial dynamics are dominant. In a perching bat, a transition from one
mode to another triggered by propreoceptive sensing and neural reflexes occur
at a remarkable high speed.

This hybrid model will allow the individual examination of each mode. Fluid-
structure interactions are hard to model, however, the aerial flip turns
(backspins) are mainly the result of differential inertial forces and are
rigorously explainable by mathematical models. To describe mode (1), the
current schemes can take a number of experimentally validated and successful
forms [3]. Therefore, the focus will be mode (2) and that what would be the
best strategy to prepare for landing on steel frames upside-down.

One may raise the concern that how this transition – i.e., mode (1) to (2) –
will occur in the first place, or, how the aerodynamic forces could vanish
quickly? I will draw inspiration from biological examples and that there are a
variety of deployable morpho-functional systems to consider [4]. In bats, the
transition occur similar to birds by facing the inside of the wings towards
incoming flow to initiate stalled conditions at which point the lift force
vanishes followed by swiftly collapsing the wings at the onset of the landing
maneuver, which almost likely completely zeros the effects of delayed stall
and aerodynamic resistance during the backspin motion.

Without such a mediolateral movement in wings, at high angle of attack
maneuvers the aerodynamic forces will be roughly normal to the wing surface at
all times resulting in large pressure forces that dominate the shear viscose
forces acting parallel to the wings. Because fixed-wing systems have no
retraction mechanism in their aerodynamic control surfaces they cannot easily
backspin. Rotory-wing systems to be able to preform backspin, they require
thruster units that can generate positive and negative aerodynamic force
relative to the body and that can be achieved through a number of costly
approaches. Another issue concerning fluid-structure interaction can arise
here. The ”ceiling effects” are unknown for us but I believe this phenomena
can cause instability in a similar way that ground effects can lead to control
issues and instability for near-ground hovering rotory-wing systems.

An informal comparison of a landing bat and a perching rotorcraft is
meaningfully relatable to the comparison between the performance of a pitcher
and a robotic arm in throwing the baseball in a baseball game. One employs
natural dynamics to perfection to pitch a ball the other is carefully
restricted to its predefined joint trajectories completely suppressing the
natural dynamics. Of course, unlike the fully-actuated manipulator, a
quadrotor is underactuated and its internal dynamics is invariant of the
supervisory controller which allows for a limited contribution from the
natural dynamics.

##  III Brief Overview of ”Harpoon” System

The hypothesis explained earlier is tested with the help of a self-sustained
ballistic object called ”Harpoon” (shown in Fig. 1). This robot is considered
as the landing gear in a morpho-functional machine. Here, the objective is to
reorient the robot towards an imaginary ceiling surface and extend (or shoot)
the landing-gear towards the ceiling. This scenario is designed after
carefully observing the high-speed images of bat landing maneuvers and once
successfully reconstructed it will capture the dynamics once inspected by [2,
5, 6, 7].

Harpoon platform is laser-cut from a thin carbon fiber plate and hosts: two
actuators, two motor controller, a receiver and a processor. The overall
design allows for a limited self-sustained sensing, actuation, computation and
communication between the robot and master computer. A brush-less motor
adjusts the angle of a bob (appendage) with respect to the body. This bob is
designed to precisely capture the body to appendage ratio found in bats. A
trigger mechanism releases the landing gear once it is faced towards the
surface. Since the projectile is charged with massive elastic energy beyond
the output performance of the trigger actuator, a mechanism with 1:40
mechanical advantage is designed and fabricated at sub-millimeter-scale
resolution as shown in Fig. 1.

The orientation angles including the Euler parameters roll, pitch and yaw are
estimated at 200 Hz using an OptiTrack system with six cameras in a 8 m3 space
and are fed back to a discrete controller with known asymptotic stability
properties. Then, the computed control inputs are fed back to the robot in a
wireless communication bus in Ethernet frames. This architecture allows for a
robust communication for fast sensing and actuation. The computations take
place off-line and it allows for an expensive real-time processing otherwise
impossible with the limited on-board processing capabilities.

Additionally, a computer-automated mechanism, called releaser, ensures the
robot is released with zero-angular momentum, therefore, the ballistic motion
will be of zero-angular-momentum nature [8]. This condition may put
restrictions on the boundaries of the trajectories from mode (1), however, a
change in the mechanism design can easily provide non-zero-angular-momentum
turns, which can easily allow for the study of non-zero-angular-momentum
backspins. To completely suppress the aerodynamic effects, wings are excluded
from the model and the host morpho-functional machine will be unveiled in the
future works.

##  IV Manipulation of inertial dynamics

The system dynamics and the equations of motion are described in a coordinate-
free fashion, i.e., a system on manifold without considering any local
coordinate charts. The kinetic energy is considered as a Riemannian metric and
in writing the discrete Hamilton’s principle associated Riemannian inter-
connections are considered. The discrete-time model on manifold with body-
fixed forces and invariant kinetic energies is erected based on the Lie group
variational integrator evolving on $SE(3)\times SO(3)$ which unlike a rigid
body results in a nontrivial dynamical model.

The numerical integrator obtained from the discrete variational principle
exhibits excellent geometric conservation properties [9] and because of this
it is used here. This means that the robot attitude automatically evolves on
the rotation group embedded in the space of special orthogonal matrices
whereas the required angular velocity of appendage is computed at the level of
the Lie algebra and the matrix exponential are employed to update the
reference solutions.

Last, an almost-globally stabilizing, geometric discrete-time controller based
on [10] will be applied to regulate the attitude in preparation for landing.

###  IV-A Notation and Assumptions

Figure 2: A hybrid model is considered to explain the landing maneuver and
Moser-Veselov description of the discrete-time inertial manipulation is
considered on $SE(3)\times SO(3)$.

It is assumed that a body coordinate frame is fixated to the robot and it
position, orientation and angular velocity are denoted by $p$, $R_{b}$ and
$\Omega_{b}$, respectively. Another body coordinate frame is attached to the
appendage where the orientation and angular velocities are denoted by $R_{a}$
and $\Omega_{a}$, respectively. Therefore, the configuration space (C-space)
is $\left\\{p,R_{b},R_{a},\Omega_{b},\Omega_{a}\right\\}$. The physical
properties including the inertia matrix and mass are denoted by $J_{a,b}$ and
$m_{a,b}$. A small change in $R_{i}$ is denoted by $F_{i}$ and $\eta_{i}$
denotes variations on $SO(3)$. The Lagrangian functional, kinetic and
potential energy are denoted by $L$, $K$ and $V$, respectively. Other widely
used notations such as $\hat{(.)}$ for wedge operator, $X[k]$ for a discrete
value of $X$ at the kth sample, $\mathrm{tr}(.)$ for the trace of a matrix,
$(.)^{T}$ as the transpose operator, $\mathrm{eig}(.)$ as the eigen-values of
a matrix and $\delta(.)$ as the variation operator are adopted. The control
input is denoted by $u$, the identity matrix is denoted by $I$, $e_{3}$ is the
z-axis unit vector and $r_{a,b}$ are the distance between the joint and the
center of mass in body and appendage.

###  IV-B Moser-Veselov Description of the Discrete-time, Constrained Dynamics

Sensing and actuation occur in discrete-time, therefore, the discrete-time
version of the dynamical system on $SE(3)\times SO(3)$ is considered and that
arises the constrained dynamics once investigated by [11]. Indirect methods
based on Rodrigues’ approach to approximate matrix exponential and first-order
approximations of variations on Lie groups are considered as well which will
be explained in this section.

####  IV-B1 Schur-form update law to $R_{b}[k]$ and $R_{a}[k]$

The discretized version of Hamilton’s principle determines the equations of
motion from a variational problem. The Lie group variational dynamics for the
system of body and appendage when external forces are absent are given by:

$\displaystyle h\hat{\Pi}_{b}[k]$
$\displaystyle=F_{b}[k]J_{b}-J_{b}F_{b}^{T}[k],$ (1) $\displaystyle
R_{b}[k+1]$ $\displaystyle=R_{b}[k]F_{b}[k],$ (2) $\displaystyle\Pi_{b}[k+1]$
$\displaystyle=F_{b}^{T}[k]\Pi_{b}[k]+hu,$ (3) $\displaystyle 0$
$\displaystyle\leq\left(h\hat{\Pi}_{b}[k]\right)^{2}+4\left(\frac{1}{2}\mathrm{tr}(J_{b})I-J_{b}\right)^{2}$
(4)

where $u\in\mathbb{R}^{3}$ allows for the inertial dynamics contributions from
the appendages. In Eqs. 1, 2, 3 and 4, the solution to the discrete-time
dynamical system can be resolved first by obtaining an answer to $F_{b}[k]$,
which involves solving the continuous-time algebraic Riccati equation (ARE)
given by Eq. 1. However, the solutions to the ARE are not unique and for these
solution to belong to special orthogonal solutions ($F_{b}[k]\in SO(3)$) the
convex constraint given by Eq. 4 must be satisfied [11]. In general, a change
in the orientation of the body $F_{b}[k]$ cannot be uniquely mapped to
$\Omega_{b}$ as $F_{b}[k]$ are minuscule changes in $R_{b}[k]\in SO(3)$. To
obtain $F_{b}[k]$, the steps from [12] are applied. The Hamiltonian matrix is
given by

$H[k]=\left[\begin{array}[]{cc}\frac{h\hat{\Pi}_{b}[k]}{2}&I\\\
\frac{\left(h\hat{\Pi}_{b}[k]\right)^{2}}{4}+J_{b}^{2}&\frac{h\hat{\Pi}_{b}[k]}{2}\end{array}\right]$
(5)

and after synthesizing the Schur form of the Hamiltonian Matrix
$\gamma[k]^{T}H[k]\gamma[k]=\left[\begin{array}[]{cc}\alpha_{1}&\alpha_{2}\\\
0&\alpha_{3}\end{array}\right]_{\mathrm{Schur},k}$, $\gamma[k]$ is partitioned
as $\gamma[k]=\left[\begin{array}[]{cc}\gamma_{1}&\gamma_{2}\\\
\gamma_{3}&\gamma_{4}\end{array}\right]_{k}$. Hence, the special orthogonal
solutions to the ARE are given by

$F_{b}[k]=\left(\gamma_{3}\gamma_{1}^{-1}+\frac{h\hat{\Pi}_{b}[k]}{2}\right)J_{b}^{-1}.$
(6)

####  IV-B2 Rodrigues’ approximation to the variations in $R_{b}[k]$ and
$R_{a}[k]$

In addition to the Schur-form solutions, another approach is considered. In
this way, special orthogonal solutions to Eqs. 1-4 are iteratively resolved by
expressing $F_{i}[k]\in SO(3),~{}~{}i\in{a,b}$ as an exponential form of
$\hat{f}_{i}[k]\in\mathfrak{so}(3)$ where $f_{i}[k]$ belongs to
$\mathbb{R}^{3}$. The Rodrigues’ formula for the changes in the body
orientation gives

$\displaystyle F_{b}[k]$
$\displaystyle=I+\frac{\sin\|f_{b}[k]\|}{\|f_{b}[k]\|}\hat{f}_{b}[k]$
$\displaystyle+\frac{1-\cos\|f_{b}[k]\|}{\|f_{b}[k]\|^{2}}\left(\hat{f}_{b}[k]\right)^{2}$

and the ARE given in Eq. 1 is re-written into the equivalent vector form given
by

$\displaystyle h\Pi_{b}[k]$
$\displaystyle=\frac{\sin\|f_{b}[k]\|}{\|f_{b}[k]\|}J_{b}f_{b}[k]$
$\displaystyle+\frac{1-\cos\|f_{b}[k]\|}{\|f_{b}[k]\|^{2}}\hat{f}_{b}[k]J_{b}f_{b}[k]$

Rodrigues’ approach was considered in addition to the Schur method to
provision for the numeric difficulties endured when attempted to use Eqs. 1-4
in a real-time reference trajectory governor. The Eqs. 1-4 yield the
symplectic geometry of the problem in hand and is followed by computational
issues about solutions in finding invariant subspaces of a Hamiltonian
matrics. Other methods that could take into account the particular structure
of Hamiltonian matrices are not less costly. Worthy of noting is that it was
noticed – in this experimentally motivated work – extra care must be paid to
as not all of the sampling rates (time intervals denoted by $h$) can guarantee
the existence and uniqueness of the solutions in Eqs. 1-4. Intuitively, this
makes sense as a relatively large time interval leads to an approximations of
$F_{i}[k]$ that violates the topological structure of $SO(3)$. The following
equation can help observe this problem

$\displaystyle F_{i}[k]=h\hat{\Omega}_{i}[k]+I$ (7)

Obviously, $h$ cannot take any arbitrary number. This equation is the direct
result of defining $F_{i}[k]$ such that $R_{i}[k+1]=R_{i}[k]F_{i}[k]$. Using
the kinematic relationship and considering an approximation for
$\hat{\Omega}_{i}[k]\approx
R^{T}_{i}[k]\left(\frac{R_{i}[k+1]-R_{i}[k]}{h}\right)$ Eq. 7 is obtained.

Almost all of these methods resolve unique positive definite solutions for ARE
when sample intervals are small, however, showing how the size of the sample
intervals can violate the Hamiltonian matrix with no imaginary eigenvalues is
well beyond the scope of this work.

Although the sufficient and necessary convex quadratic constraint in Eq. 4 is
a strong condition that secures the existence of such a solution, its
dependence to the sampling time intervals can be limiting.

When $\left(h\hat{\Pi}_{b}[k]\right)^{2}/4+J_{b}^{2}$ is positive definite
then Eq. 1 has a unique positive semi-definite solution $S[k]\geq 0$, i.e.,
$\mathrm{eig}(S[k]+\hat{\Pi}_{b}[k]/2)$ has positive real parts. Consequently,
$F_{b}[k]=(S[k]+h\hat{\Pi}_{b}[k]/2)J_{b}^{-1}$ is the special orthogonal
matrix being sought. When $\left(h\hat{\Pi}_{b}[k]\right)^{2}/4+J_{b}^{2}$ is
positive definite an invertible matrix $E[k]$ exists such that
$\left(h\hat{\Pi}_{b}[k]\right)^{2}/4+J_{b}^{2}=E[k]^{T}E[k]$. Consequently,
$(h\hat{\Pi}_{b}[k]/2,I)$ and $(E[k],h\hat{\Pi}_{b}[k]/2)$ will be
stabilizable and detectable, respectively, which results in
$\mathrm{eig}(S[k]+h\hat{\Pi}_{b}[k]/2)$ with positive real parts. Since
$\mathrm{eig}(\pm h\hat{\Pi}_{b}[k]/2-I)$ are the form $-1\pm\beta_{i}$ they
have negative real part. Therefore, $(h\hat{\Pi}_{b}[k]/2+I)$ and
$(-h\hat{\Pi}_{b}[k]/2,E^{T})$ are stabilizable, that is,
$(h\hat{\Pi}_{b}[k]/2,I)$ is stabilizable and $(E,h\hat{\Pi}_{b}[k]/2)$ is
detectable.

####  IV-B3 Base-line coordinate-free system on $SO(3)$ manifold

Based on Hamilton’s Principle, the action integral is given from the
Lagrangian $L=K-V$ and its variations are zero. Because the configuration
group is the Lie group, this variation should be consistent with this
geometry, therefore, the varied $\Omega_{i},i\in{a,b}$ in continuous-time mode
is given by

$\delta\Omega_{i}=\hat{\Omega}_{i}\eta_{i}+\dot{\eta}_{i},~{}i\in{a,b}\\\ $

and in discrete-time mode is given by

$\delta\Omega_{i}[k]=\frac{1}{h}\left(F_{i}[k]\eta_{i}[k+1]-\eta_{i}[k]\right),~{}i\in{a,b}.\\\
$

Consequently the variations of the action integral are obtained. In driving
the variations of continuous-time action integral,
${\Omega_{i}^{T}}{J_{i}}{\dot{\eta}_{i}}={\frac{d}{dt}}({\Omega_{i}^{T}}{J_{i}}{\eta_{i}})-{\frac{d{\Omega_{i}^{T}}}{dt}}{J_{i}}{\eta_{i}}$
is used. Note that the continuous- and discrete-time variations vanish at the
time boundaries, that is, $t_{0}$ and $t_{f}$ for the continuous system and
$k=0$ and $k=N$ for the discrete system. From Hamilton’s principle, the above
continuous- and discrete-time variations in the action integral should be zero
for all variation of Lie algebra. Hence, the continuous time system on
manifold is given by:

$\dot{x}=G_{1}(x)+G_{2}(x)u$ (8)

where $x=(p;R_{b};R_{a};\dot{p};\Omega_{b};\Omega_{a})\in SE(3)\times SO(3)$
and $u\in\mathbb{R}^{3}$ are the states and control inputs, respectively;
$G_{1}(x)=\psi^{-1}(x)\phi(x)$ and $G_{2}(x)=\psi^{-1}(x)\left[0_{3\times
3},I_{3\times 3},-R_{a}^{T}R_{b}\right]^{T}$ are obtained after sorting the
varaitions and applying some basic algebra. The symmetric matrix
$\psi(x)\in\mathbb{R}^{9\times 9}$ is given by

$\psi(x)=\left[\begin{array}[]{ccc}-MI&\xi_{b}\nu_{b}&-\xi_{a}\nu_{a}\\\
\xi_{b}\nu^{T}_{b}&r_{b}\xi_{b}\hat{e}^{2}_{3}-J_{b}&r_{a}\xi_{b}\nu^{T}_{b}\nu_{a}\\\
-\xi_{a}\nu^{T}_{a}&r_{a}\xi_{b}\nu_{a}^{T}\nu_{b}&r_{a}\xi_{a}\hat{e}^{2}_{3}-J_{a}\end{array}\right]$
(9)

and $\phi(x)\in\mathbb{R}^{9}$ is given by

$\phi(x)=\left[\begin{array}[]{l}\xi_{b}R_{b}\kappa_{b}-\xi_{a}R_{a}\kappa_{a}-Mge_{3}\vspace{0.25cm}\\\
\hat{\Omega}_{b}\Pi_{b}+2\xi_{b}\zeta^{T}_{b}\dot{p}+r_{b}\xi_{b}\hat{e}_{3}\kappa_{b}+\dots\\\
2r_{b}\xi_{a}\zeta^{T}_{b}\nu_{a}\Omega_{a}+\xi_{b}g\hat{\nu_{b}}e_{3}-r_{b}\xi_{a}\nu_{b}\zeta_{a}\Omega_{a}\vspace{0.25cm}\\\
\hat{\Omega}_{a}\Pi_{a}-2\xi_{a}\zeta^{T}_{a}\dot{p}+r_{a}\xi_{a}\hat{e}_{3}\kappa_{b}-\dots\\\
2r_{a}\xi_{b}\zeta^{T}_{a}\nu_{b}\Omega_{b}-\xi_{b}g\hat{\nu_{a}}e_{3}-r_{a}\xi_{b}\nu_{a}\zeta_{a}\Omega_{b}\end{array}\right]$
(10)

where $M=\sum_{i}m_{i}$, $\xi_{i}=m_{a}r_{i}$, $\nu_{i}=R_{i}\hat{e}_{3}$,
$\zeta_{i}=R_{i}\hat{\Omega}_{i}\hat{e}_{3}$ and
$\kappa_{i}=\hat{\Omega}^{2}_{i}e_{3}$ for $i\in{a,b}$. The discreteization of
this continuous dynamics (Eq.1-4) is achieved by applying the discrete-time
variation explained above and applying Hamilton’s principle. Next, I will
briefly describe the closed-loop feedback used to regulate the body
orientation.

###  IV-C Discrete feedback control

A scalar measure of the attitude error is given by the rotation angle
$\theta[k]$ about the eigen-axis needed to rotate the robot form $R_{b}[k]$ to
the desired attitude $R_{d}[k]$ and is given by
$\theta[k]=\cos^{-1}\left(\frac{1}{2}\left[\mathrm{tr}R^{d}_{T}[k]R_{R}[k]-1\right]\right)$.
This metric is applied to indicate the time of launching the landing gear.

The controller applied here does not require any knowledge in regard to the
mass and inertia of the robot body and appendage. A major limitation in
designing the closed-loop feedback for this constrained dynamics is that not
all of the arbitrary reference trajectories can be enforced by the tracking
controller and the tracking problem is subject to the constraints given by Eq.
4.

The issue of constraint satisfaction in nonlinear systems affine- and non-
affine-in-control and its separation from the issue of closed-loop system
design have been addressed very extensively. Almost all of these works do so
by introducing a reference governor, an auxiliary nonlinear device that
operates between the reference command and the input to the closed-loop
system. A similar approach to [13] is adopted here.

The control input with known stability properties [10] is given by

$\displaystyle u=K_{1}\Pi_{b}[k]-e[k]$ (11)

where $K_{1}$ is the feedback gain and $e[k]$ can be calculated by the
following relation,

$\displaystyle\hat{e}[k]=K_{2}\Gamma^{T}[k]R_{b}[k]-R_{b}^{T}[k]\Gamma[k]K_{2}$
(12)

where $K_{2}$ is a symmetric positive-definite matrix and $\Gamma[k]$ is the
reference input to the body orientation and is resolved by finding $\alpha[k]$
in a nonlinear optimizer subject to the following constraints:

$\displaystyle 0$ $\displaystyle\leq\alpha[k],\vspace{0.25cm}$ (13)
$\displaystyle\alpha[k]$ $\displaystyle\leq 1,\vspace{0.25cm}$ (14)
$\displaystyle 0$
$\displaystyle\leq\left(h\hat{\Pi}_{b}[k]\right)^{2}+4\left(\frac{1}{2}\mathrm{tr}(J_{b})I-J_{b}\right)^{2},\vspace{0.25cm}$
(15) $\displaystyle 0$
$\displaystyle=\Gamma[k]-\left(R_{d}[k]\Gamma^{T}[k-1]\right)^{\alpha[k]}\Gamma[k-1].$
(16)

Note that $R_{d}[k]$ is the desired body orientation before considering the
convex constraint given by Eq. 4.

##  V Numerical & Experimental Results

The appendage inertial model is scaled to account for 16% of the total body
mass according to [7, 14]. A total estimated body weight of 30 gr is
considered, which is obtained after assuming a density of $2\times
10^{3}kg.m^{-3}$ according to [14], while equally distanced bobs from the body
joint (shown in Fig. 2) are considered. A simple body geometry is considered
for a seamless inertia matching.

When only the dynamics restricted to the longitudinal plan is considered the
C-space is the two-sphere $S^{2}=\left\\{q\in\mathbb{R}^{3}|q.q=1\right\\}$,
where $q$ is a unit vector (Note that the body and appendage angles are
embedded in $q$). The tangent space at $q$, namely $T_{q}S^{2}$, is identified
with $T_{q}S^{2}\approx\left\\{\omega\in\mathbb{R}^{3}|q.\omega=0\right\\}$
where $\omega\in\mathbb{R}^{3}$ represents the angular velocity.

A fixed desired body orientation $R_{d}$ is constructed according to the Euler
convention by ramping the pitch angles to 120 and 180 deg and updated by the
reference governor subject to the constraints explained above (shown in Fig.
4(a) and 4(b)). Samples (purple dots) are considered by the optimizer at the
neighborhood of the desired angles (black dots) before resolving the actual
trajectories. Note that due to issues in the cluttered raster image some of
the purple dots are plotted. In Fig. 5(a), the orientations of a 3D-rigid body
representing the body (note that the appendage is not shown) are shown along
the time axis for a simultaneous tracking of roll and pitch angels -20 and 60
deg, respectively. In Fig. 5(b), the simulated Euler angles are shown.

In experiment, only the longitudinal motion is considered. The framework was
tested in a series of experiments devised with an off-board computer, motion
capture system and high-speed imaging equipment. In Fig. 3, the snapshots from
one experiment are shown. The computer-automated releaser mechanism releases
the robot with zero-angular momentum in a free fall configuration while the
motion capture system captures the body and appendage orientations. The
commanded $\Omega_{a}$ (shown in Fig. 6) are resolved and transmitted to the
robot for the desired body orientation angle 180 deg.

Figure 3: Snapshots of Harpoon while reorienting towards the virtual ceiling
to launch the landing gear

(a) Resolved reference trajectory (green) reaching to 120 deg subject to the
constraints in Eq. 4 shown on a two-sphere

(b) Resolved reference trajectory reaching to 180 deg

(a) Simulated body motion

(b) Simulated appendage movements to achieve $-20$ and $60$ degrees variations
in the roll and pitch angles

Figure 6: Appendage angular velocity $\Omega_{ay}$

##  VI Conclusion

The ability of animals to reorient their bodies mid-air without the benefit of
net external torques is known yet remarkably has been overlooked in aerial
drone design. Instead, current designs consider thrust vectoring for agile
aerial robotics, which is known to be limited based on aerodynamics laws. This
work used a tiny robot called Harpoon with prohibitive restrictions associated
with size and payload not allowing for multi-thruster designs typically found
in quadrotors and demonstrated that extremely fast body reorientation and
preparation for upside-down landing unique to bats is possible through the
closed-loop manipulations of inertial dynamics in the robot. To prepare for
landing, a rubber-band-propelled landing gear design and associated trigger
mechanism was proposed. The closed-loop manipulations of inertial dynamics
took place based on a symplectic description of the dynamical system (body and
appendage), which is known to exhibit an excellent geometric conservation
properties.

## References

  * [1] D. K. Riskin, A. Bergou, K. S. Breuer, and S. M. Swartz, “Upstroke wing flexion and the inertial cost of bat flight,” _Proceedings of the Royal Society B: Biological Sciences_ , vol. 279, no. 1740, pp. 2945–2950, 2012.
  * [2] A. J. Bergou, L. Ristroph, J. Guckenheimer, I. Cohen, and Z. J. Wang, “Fruit flies modulate passive wing pitching to generate in-flight turns,” _Physical review letters_ , vol. 104, no. 14, p. 148101, 2010.
  * [3] A. Ramezani, S.-J. Chung, and S. Hutchinson, “A biomimetic robotic platform to study flight specializations of bats,” _Science Robotics_ , vol. 2, no. 3, pp. Art–No, 2017.
  * [4] J. F. V. Vincent, “Deployable Structures in Biology,” in _Morpho-functional Machines: The New Species_ , F. Hara and R. Pfeifer, Eds. Springer Japan, 2003, pp. 23–40.
  * [5] D. K. Riskin, J. W. Bahlman, T. Y. Hubel, J. M. Ratcliffe, T. H. Kunz, and S. M. Swartz, “Bats go head-under-heels: the biomechanics of landing on a ceiling.” _The Journal of experimental biology_ , vol. 212, no. Pt, pp. 945–953, 2009.
  * [6] Bahlman Joseph W., Swartz Sharon M., Riskin Daniel K., and Breuer Kenneth S., “Glide performance and aerodynamics of non-equilibrium glides in northern flying squirrels (Glaucomys sabrinus),” _Journal of The Royal Society Interface_ , vol. 10, no. 80, p. 20120794, Mar. 2013.
  * [7] J. Iriarte-Diaz and S. M. Swartz, “Kinematics of slow turn maneuvering in the fruit bat Cynopterus brachyotis.” _The Journal of experimental biology_ , vol. 211, no. Pt, pp. 3478–3489, 2008.
  * [8] A. Arabyan and D. Tsai, “A distributed control model for the air-righting of a cat,” p. 9.
  * [9] E. Hairer, C. Lubich, and G. Wanner, _Geometric numerical integration: structure-preserving algorithms for ordinary differential equations_. Springer Science & Business Media, 2006, vol. 31.
  * [10] A. Sanyal, A. Fosbury, N. Chaturvedi, and D. Bernstein, “Inertia-free spacecraft attitude tracking with disturbance rejection and almost global stabilization,” _Journal of guidance, control, and dynamics_ , vol. 32, no. 4, pp. 1167–1178, 2009.
  * [11] J. Moser and A. P. Veselov, “Discrete versions of some classical integrable systems and factorization of matrix polynomials,” _Communications in Mathematical Physics_ , vol. 139, no. 2, pp. 217–243, 1991.
  * [12] A. Bunse-Gerstner and V. Mehrmann, “A symplectic qr like algorithm for the solution of the real algebraic riccati equation,” _IEEE Transactions on Automatic Control_ , vol. 31, no. 12, pp. 1104–1113, 1986.
  * [13] E. Gilbert and I. Kolmanovsky, “Nonlinear tracking control in the presence of state and control constraints: a generalized reference governor,” _Automatica_ , vol. 38, no. 12, pp. 2063–2073, 2002.
  * [14] J. Iriarte-Diaz, D. K. Riskin, D. J. Willis, K. S. Breuer, and S. M. Swartz, “Whole-body kinematics of a fruit bat reveal the influence of wing inertia on body accelerations,” _Journal of Experimental Biology_ , vol. 214, no. 9, pp. 1546–1553, May 2011.