# Machine Learning Regression based Single Event Transient Modeling Method for
Circuit-Level Simulation

ChangQing Xu  Member, IEEE  Yi Liu  XinFang Liao  JiaLiang Cheng and YinTang
Yang  Member, IEEE  “This work was supported in part by the Fundamental
Research Funds for the Central Universities, by the Industry-University-
Academy Cooperation Program of Xidian University-Chongqing IC Innovation
Research Institute under Grant CQIRI-2021CXY-Z01, by Equipment Pre-research
Project of China under Grant 41424050607.” Changqing Xu,Yi Liu, XinFang Liao,
JiaLiang Cheng, YinTang Yangon, are with the School of Microelectronics,
Xidian University, Xi’an 710071, China (e-mail: cqxu@xidian.edu.cn).

###### Abstract

In this paper, a novel machine learning regression based single event
transient (SET) modeling method is proposed. The proposed method can obtain a
reasonable and accurate model without considering the complex physical
mechanism. We got plenty of SET current data of SMIC 130nm bulk CMOS by TCAD
simulation under different conditions (e.g. different LET and different drain
bias voltage). A multilayer feedfordward neural network is used to build the
SET pulse current model by learning the data from TCAD simulation. The
proposed model is validated with the simulation results from TCAD simulation.
The trained SET pulse current model is implemented as a Verilog-A current
source in the Cadence Spectre circuit simulator and an inverter with five fan-
outs is used to show the practicability and reasonableness of the proposed SET
pulse current model for circuit-level single-event effect (SEE) simulation.

linear energy transfer(LET), machine learning, neural network, single event
transient(SET).

##  1 Introduction

With continuous downscaling of transistor sizes and supply voltages, the
sensitivity to single event transient (SET) has become one of the most
important reliability issues for aerospace integrated circuit (IC) which is
under radiation exposure [1, 2, 3]. Generally, the evaluation and analysis of
the circuit’s sensitivity to SETs can be performed either by the computer-
aided modeling and simulations or through the experimental approach. As the
irradiation experiments are very expensive, circuit-level single event
transient (SET) simulation is suitable for understanding how the circuit
responds to the ion strike and how the circuit is hardened to avoid SET at the
early stages of the design [6]. SET model has been widely applied in analogy
and digital circuit SET simulation[4, 5, 6], which can guide radiation-
hardened IC design.

Figure 1: Flow of the SET modeling method based on Machine Learning
Regression.

However, the accuracy of the SET model have a direct impact on circuit-level
SET simulation.To obtain an accurate SET model, lots of researchers proposed
various approaches to build SET model [7, 8, 9, 10, 11]. A double-exponential
model [2, 7] is widely used to model the SET current from the drain to source
of NMOS. Another exponential model is proposed by Freeman [8], in which a
single technology-related timing parameter is defined to make the mode
suitable for different technologies. Hu’s current model considering the angel
of particle incidence [9]. Considering the influence of strike location, a
bias-dependent single event model is proposed in which a drift factor and a
diffusion factor is introduced to produce the accurate single-event current
pulse [11]. There are also some physics-based SET models [12, 13], which can
transient potential and drain current on linear energy transfers (LETs),
strike positions, drain and gate biases, etc. Physics-based SET can capture
the essential physics behind these current transients, but these models are
too complex to use in circuit-level simulation. Those SET models, which are
based on certain functions, are suitable for circuit-level simulation.
However, due to the complex physical mechanism of devices and single event
effect, it is hard to build an accurate SET current model based on a certain
function, such as double-exponential function, when various variables, such as
LET, strike location, bias voltage, etc, are considered.

In this paper, we proposed a novel SET modeling method based on machine
learning regression technology. A feedfordward neural network is used to train
the SET pulse current model. The trained model is validated and implemented as
a Verilog-A current source in the Cadence Spectre circuit simulator.

##  2 Machine Learning Regression based SET Modeling Method

Fig. 1 shows the flow of SET modeling method based on machine learning
regression where all steps have been partition into phases. The Cogenda
VisualTCAD[14], which consists of Gds2mesh, VisualParticle, and genius, is
used to generate a dataset for training and testing the SET pulse current
model. Simulation results are collected over various LET, drain bias voltage,
and simulation time. These data will be pre-processed and partitioned as
training and testing datasets. A multi-layer feedforward neural network is
applied as SET model and trained based on the dataset obtained from TCAD
simulation. To make a trade-off between the accuracy and the complexity of the
SET model, the parameters, such as neural network layers, number of neurons,
transfer functions, etc, are tuned based on the accuracy requirement.

###  2.1 SET Simulation in TCAD and Data Collection

Figure 2: 3D structure of the NMOS.

In this paper, all TCAD simulation results are obtained by Cogenda
VisualTCAD[14]. Fig. 2 shows the 3D structure of the NMOS in a Semiconductor
Manufacturing International Corporation (SMIC) 130nm CMOS technology. The
device dimension of NMOS for the model and TCAD simualtions in the following
sections are channel length $L$ = 130 nm, channel width ${W}$ = 420 nm, gate
oxide thinkness ${t_{ox}}$ = 2.58 nm, source/darin doping is $10^{20}cm^{-3}$,
and channel doping $N=1.57\times 10^{17}cm^{-3}$. The gate voltage is fixed at
0 V, source and substrate are grounded and the drain voltage is set to $0\sim
1.8V$.

The heavy ion trike is simualted using heavy ion model in Cogenda VisualTCAD
[14]. The models used in the simulation are Shockley-Read-Hall, Auger, and
direct (or radiative) recombination. SRH carrier lifetime has both doping
dependence and temperature dependence. Other models applied are effective
density of states model, Fermi–Dirac statistics, impact ionization, and Lucent
mobility model. Heavy ion samples with LET of $4\sim 100MeV\cdot cm^{2}/mg$,
The incident positions of these ions were the NMOS drain center, the track
length of the heavy ions was set to be 10 $\mu m$, and the incident direction
is the direction vertical to the surface of the device. TCAD simulation was
conducted for SMIC 130nm NMOS and the drain current data at various drain bias
voltages($0\sim 1.8V$) and various LETs($4\sim 100MeV\cdot cm^{2}/mg$) is
collected. LET, drain bias voltage and time have been considered as input
features and drain current has been divided as output features. Because the
data obtained in TCAD simulation is too dense at the corner of the curve which
may cause neural network over fitting, cubic spline data interpolation is used
to get the sampling point based on accuracy requirements. 418000 data is
obtained based on TCAD simulation results. Among these data, 70% of data will
be selected randomly and considered as training data, 15% data will divided
into validation and remaining 15% data will be set as test data.

###  2.2 SET pulse current model based on feedfordword neural network

In this paper, feedforward neural networks, which consist of an input layer,
several hidden layers, and an output layer, are used to build SET pulse
current model. Levenberg–Marquardt algorithm is used as learning rule. Because
the number of neurons, number of layers and transfer function applied will
directly affect the fitting ability of the neural network. Several neural
networks are trained 1000 epoch to build SET pulse current model and the mean
square errors(MSE) are shown in Table 1. Considering the complexity and
accuracy of the SET model, the neural network(8$\times$8$\times$1) is selected
to build SET pulse current model.

Table 1: MSE of different neural networks after training Network | Tranfer Function | Error (MSE)
---|---|---
Hidden Layer | Output Layer
16$\times$1 | tansig | purelin | 9.54e-3
16$\times$1 | logsig | purelin | 2.48e-2
16$\times$1 | elliotsig | purelin | 1.83e-2
8$\times$8$\times$1 | tansig | purelin | 1.37e-3
8$\times$8$\times$1 | logsig | purelin | 2.14e-3
8$\times$8$\times$1 | elliotsig | purelin | 2.18e-2
8$\times$16$\times$8$\times$1 | tansig | purelin | 2.76e-4
8$\times$16$\times$8$\times$1 | logsig | purelin | 4.99e-4
8$\times$16$\times$8$\times$1 | elliotsig | purelin | 2.10e-2

Figure 3: SET model based on feedfordward neural network

Fig. 3 shows SET pulse current model based on 8$\times$8$\times$1 feedfordward
neural network. In the input layer, time, LET and bias voltage are the input
variables. In the hidden layers, tansig function is used as transfer function.
In the output layer, a simple line function is selected as transfer function.
$W^{(1)}_{8\times 3}$, $W^{(2)}_{8\times 8}$, $W^{(3)}_{1\times 8}$ are
trained weight matrices for corresponding layers. $b^{(1)}_{8\times 1}$,
$b^{(2)}_{8\times 1}$, $b^{(3)}_{1\times 1}$ are trained bias matrices for
corresponding layers. Due to different input data units, all input data need
to be normalized first and the neural network output will be denormalized to
generate the accurate SET pulse current.

##  3 Results and Discussion

Figure 4: SET current versus time for different drain voltage and LET

Figure 5: SET current and output voltage of the two-stage inverter versus time
for different LET

To show the accuracy of the trained SET pulse current model, we compare
trained model and TCAD simulation results. Fig. 4 shows the SET pulse current
model under various drain bias voltage (0.2, 0.4, 0.6, 0.8, 1.2V) and various
LET (5, 20, 40, 60, 80 $MeV\cdot cm^{2}/mg$) with respect to time and compared
with TCAD simulation. As Fig. 4 is shown, the trained model can fit the SET
pulse current very well, especially for the peak value and width of current
pulse which is of great importance to circuit-level simulation.

The developed SET pulse current model is implemented as a Verilog-A current
source $I_{SET}$ in the Cadence Spectre circuit simulator. The SET simulation
is carried out in the two-stage CMOS inverter by connecting the SET current
source in parallel to the OFF-state NMOS in the first-stage of inverter as
shows in Fig. 5. The first-stage inverter, which is five fin-outs, is
connected to five inverters. The ion strike at the OFF-state NMOS in the
first-stage inverter and ion strike time is 200ps. As shown in Fig. 5, the ion
strike at the OFF-state NMOS in the first-stage inverter at 200 ps, which
causes the output voltage is pulled down from the logic high state. A plateau,
which is caused by the dynamic interaction of the depressed drain voltage due
to the strike and the compensating PMOS transistor drive current [15, 16], can
be observed in the SET current of nMOS in the inverter. For higher LETs, the
voltage perturbation is more obvious due to higher amount of charge deposition
in the silicon body.

##  4 Conclusion

A novel machine learning regression based SET modeling method is proposed and
a SET pulse current model based on feedforward neural network is developed and
validated. The developed model is easy to integrate into the circuit simulator
and takes lesser computation time compared to the TCAD simulations. The
analysis can be extended to other digital or analog circuits using the
proposed SET pulse current model, which helps the circuit designers to predict
and mitigate the influence of SET on circuit design.

## References

  * [1] M. Andjelkovic, A. Ilic, Z. Stamenkovic, M. Krstic, and R. Kraemer, “An overview of the modeling and simulation of the single event transients at the circuit level,” in _2017 IEEE 30th International Conference on Microelectronics (MIEL)_ , 2017, pp. 35–44.
  * [2] L. Artola, M. Gaillardin, G. Hubert, M. Raine, and P. Paillet, “Modeling single event transients in advanced devices and ics,” _IEEE Transactions on Nuclear Science_ , vol. 62, no. 4, pp. 1528–1539, 2015.
  * [3] V. Ferlet-Cavrois, L. W. Massengill, and P. Gouker, “Single event transients in digital cmos—a review,” _IEEE Transactions on Nuclear Science_ , vol. 60, no. 3, pp. 1767–1790, 2013.
  * [4] J. S. Kauppila, L. W. Massengill, D. R. Ball, M. L. Alles, R. D. Schrimpf, T. D. Loveless, J. A. Maharrey, R. C. Quinn, and J. D. Rowe, “Geometry-aware single-event enabled compact models for sub-50 nm partially depleted silicon-on-insulator technologies,” _IEEE Transactions on Nuclear Science_ , vol. 62, no. 4, pp. 1589–1598, 2015.
  * [5] K. Foley, N. Seifert, J. B. Velamala, W. G. Bennett, and S. Gupta, “Irt: A modeling system for single event upset analysis that captures charge sharing effects,” in _2014 IEEE International Reliability Physics Symposium_ , 2014, pp. 5F.1.1–5F.1.9.
  * [6] C. Q. Xu, Y. Liu, X. D. Weng, Z. B. Li, and Y. T. Yang, “Single event transient pulses fault injection model based on let for circuit-level simulation,” in _2020 IEEE 15th International Conference on Solid-State Integrated Circuit Technology (ICSICT)_ , 2020, pp. 1–3.
  * [7] G. I. Wirth, M. G. Vieira, and F. Kastensmidt, “Accurate and computer efficient modelling of single event transients in cmos circuits,” _Iet Circuits Devices Systems_ , vol. 1, no. 2, pp. 137–142, 2007.
  * [8] L. B. Freeman, “Critical charge calculations for a bipolar sram array,” _IBM Journal of Research and Development_ , vol. 40, no. 1, pp. 119–129, 1996\.
  * [9] Hu and C., “Alpha-particle-induced field and enhanced collection of carriers,” _IEEE Electron Device Letters_ , vol. 3, no. 2, pp. 31–34, 1982\.
  * [10] T. Merelle, H. Chabane, J. . Palau, K. Castellani-Coulie, F. Wrobel, F. Saigne, B. Sagnes, J. Boch, J. R. Vaille, G. Gasiot, P. Roche, M. . Palau, and T. Carriere, “Criterion for seu occurrence in sram deduced from circuit and device simulations in case of neutron-induced ser,” _IEEE Transactions on Nuclear Science_ , vol. 52, no. 4, pp. 1148–1155, 2005\.
  * [11] L. Ding, W. Chen, T. Wang, R. Chen, Y. Luo, F. Zhang, X. Pan, H. Sun, and L. Chen, “Modeling the dependence of single-event transients on strike location for circuit-level simulation,” _IEEE Transactions on Nuclear Science_ , vol. 66, no. 6, pp. 866–874, 2019.
  * [12] Y. M. Aneesh, S. R. Sriram, K. R. Pasupathy, and B. Bindu, “An analytical model of single-event transients in double-gate mosfet for circuit simulation,” _IEEE Transactions on Electron Devices_ , vol. 66, no. 9, pp. 3710–3717, 2019.
  * [13] Y. M. Aneesh and B. Bindu, “A physics-based single event transient pulse width model for cmos vlsi circuits,” _IEEE Transactions on Device and Materials Reliability_ , vol. 20, no. 4, pp. 723–730, 2020.
  * [14] Cogenda Visual TCAD 1.9.2. [Online]. Available: {http://www.cogenda.com/article/download}
  * [15] J. S. Kauppila, A. L. Sternberg, M. L. Alles, A. M. Francis, J. Holmes, O. A. Amusan, and L. W. Massengill, “A bias-dependent single-event compact model implemented into bsim4 and a 90 nm cmos process design kit,” _IEEE Transactions on Nuclear Science_ , vol. 56, no. 6, pp. 3152–3157, 2009.
  * [16] S. DasGupta, A. F. Witulski, B. L. Bhuva, M. L. Alles, R. A. Reed, O. A. Amusan, J. R. Ahlbin, R. D. Schrimpf, and L. W. Massengill, “Effect of well and substrate potential modulation on single event pulse shape in deep submicron cmos,” _IEEE Transactions on Nuclear Science_ , vol. 54, no. 6, pp. 2407–2412, 2007.


# Dynamic and Distributed Optimization for the Allocation of Aerial Swarm
Vehicles

Jason Hughes, Dominic Larkin, Charles O’Donnell, Christopher Korpela The
authors are with the Robotics Research Center, Department of Electrical
Engineering and Computer Science, United States Military Academy, West Point,
NY, 10996 USA. E-mail: {jason.hughes, dominic.larkin, charles.o’donnell,
<EMAIL_ADDRESS>
###### Abstract

Optimal transport (OT) is a framework that can guide the design of efficient
resource allocation strategies in a network of multiple sources and targets.
This paper applies discrete OT to a swarm of UAVs in a novel way to achieve
appropriate task allocation and execution. Drone swarm deployments already
operate in multiple domains where sensors are used to gain knowledge of an
environment [1]. Use cases such as, chemical and radiation detection, and
thermal and RGB imaging create a specific need for an algorithm that considers
parameters on both the UAV and waypoint side and allows for updating the
matching scheme as the swarm gains information from the environment.
Additionally, the need for a centralized planner can be removed by using a
distributed algorithm that can dynamically update based on changes in the
swarm network or parameters. To this end, we develop a dynamic and distributed
OT algorithm that matches a UAV to the optimal waypoint based on one parameter
at the UAV and another parameter at the waypoint. We show the convergence and
allocation of the algorithm through a case study and test the algorithm’s
effectiveness against a greedy assignment algorithm in simulation.

###### Index Terms:

Discrete Optimal Transport, Resource Matching, Swarming

##  I Introduction

Optimal transport (OT) is a centralized framework that is often leveraged for
resource allocation from a set of source nodes to a set of target nodes [2].
OT has been used for many efficient resource allocation and matching problems
such as allocating raw materials for consumption, matching employees to tasks
in a corporate environment, and allocating limited power units in areas struck
by natural disasters.

While optimal transport is used for many problems, it has yet to be adapted
for mapping members of a UAV swarm to a series of waypoints in a dynamic and
distributed manner. To this end, the main contribution of this paper is a
dynamic and distributed optimal transport algorithm for the efficient mapping
of UAVs to waypoints. We consider this research as an initial development that
will lead to more complex matching problems with a heterogeneous swarm.

Leveraging the discrete optimal transport formulation has many advantages.
First, consider a network where both the UAV and the waypoint each have a
parameter to be accounted for in the optimization algorithm, as shown in Fig.
1. In this network, the distance between the UAV and the waypoint is the agent
side parameter, and the importance of visiting a waypoint is the waypoint side
parameter. With the OT algorithm, these parameters exist at the connection
between the nodes rather than at the nodes. Secondly, a dynamic algorithm is
formed to update the matching scheme when parameters in the network change.
The network parameters only change when a UAV visits a waypoint. This
constraint avoids the possibility of visiting the same waypoint twice and
prevents the algorithm from becoming stuck in a local minimum. The algorithm
also must account for when a UAV has to land for a battery swap. Finally,
consider an algorithm for application to a task where information is gained
only at a waypoint. Examples of these tasks include chemical sensing,
radiation sensing, or surveying an area.

A problem to be considered is as the network becomes larger and larger with
more nodes, the computational complexity grows exponentially. As illustrated
in Fig 1, there is a relatively small number of waypoints (40) and four UAVs,
but there are 120 connections resulting in 240 parameters. This problem
becomes more profound when the allocation occurs sequentially over time, as it
does for our matching problem. In this case, parameters are added at each
iteration, adding complexity. Since centralized control of the swarm is
undesirable, we developed a distributed algorithm using the Alternating
Direction Method of Multipliers (ADMM). This algorithm allows each UAV to
solve its own optimization problem and then communicate its results with the
other swarm members.

Figure 1: A Network $\mathcal{G}$ with four UAV nodes and forty waypoint
nodes, and a full connection scheme where every UAV is connected to every
waypoint node.

Another consideration is that the network changes as time progresses. This
change could lead to the undesirable behavior of swarm members visiting the
same waypoint more than once. There is also the problem of battery life. At
some point, UAVs may need to drop out of the swarm for a battery swap or a
malfunction. Lastly, as a UAV visits a waypoint and takes a sensor reading,
information is gained about the waypoint side of the network. The network
parameters change with this additional information, and a new optimal matching
now exists. These requirements lead to developing a dynamic distributed
algorithm that calculates an optimal solution any time the parameters and
network structure update.

Furthermore, consider the information that needs exchanging at each step. The
UAVs communicate by sending Wi-Fi packets when the algorithm runs on hardware.
This creates a need to keep the amount of data being transferred as small as
possible to support the scalability of the swarm. Once a waypoint is reached,
the proposed algorithm only needs to share three integers with the other
agents in the swarm, thus limiting the amount of data transferred.

The contributions of this paper are summarized as follows. First, a
centralized optimal transport algorithm is leveraged for the UAV to waypoints
matching problem. Second, a distributed algorithm is formed using ADMM, so
there is no need for a centralized control station. Third, a dynamic
distributed algorithm is created that can automatically calculate the optimal
solution to the matching problem when parameters are updated. Lastly, we
corroborate our results with a preliminary study in MATLAB and show the
algorithm’s effectiveness in simulation.

Related Work. Swarming UAVs have been widely studied [3], [4], [5], [6]. The
centralized and distributed optimal transport algorithm that we leverage is
developed in [7]. A dynamic and distributed optimal transport algorithm is
shown in [8], and a secure distributed OT algorithm is developed in [8].
Optimal transport has been used for swarm guidance in [9], [10] but has yet to
be used dynamically. A dynamic OT algorithm is developed in [11] but it is not
distributed. Radioactivity sensing swarms refer to waypoint allocation but do
not use an OT based algorithm [12],[13]. Distributed resource matching and
allocation algorithms are studied extensively in [14], [15], [16] and
specifically optimal transport in [2].

##  II Algorithm Formulation

In this section we set up the discrete optimal transport problem with the UAV
to waypoint matching problem.

###  II-A Discrete Optimal Transport

Start by considering each UAV as a ”target” node represented by
$x\in\mathcal{X}$ and each waypoint as a ”source” represented by
$y\in\mathcal{Y}$ and define $\mathcal{N}=\\{\mathcal{X}+\mathcal{Y}\\}$ as
the set of all nodes. Each source node is connected to some or all target
nodes, the set of targets that are connected to source $y$ is denoted by
$\mathcal{X}_{y}$, alternatively the set of source nodes connected to target
$x$ is denoted by $\mathcal{Y}_{x}$. For convenience, we denote the set of
edges connecting target and source nodes as
$\mathcal{E}:=\\{\\{x,y\\}|x\in\mathcal{X}_{y},y\in\mathcal{Y}_{x}\\}$. We
denote the network made up of the nodes and edges as
$\mathcal{G}=\\{\mathcal{N},\mathcal{E}\\}$. With this notation, we use the
optimal transport formulation from [7] that considers utility from both the
target and source sides. We have the following discrete optimal transport
formulation:

$\displaystyle\max_{\Pi}\ \sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}$
$\displaystyle
d_{xy}(\pi_{xy})+\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}s_{xy}(\pi_{xy})$
(1) $\displaystyle\mathrm{s.t.}$
$\displaystyle\underline{p}_{x}\leq\sum_{y\in\mathcal{Y}_{x}}\pi_{xy}\leq\bar{p}_{x},\
\forall x\in\mathcal{X},$
$\displaystyle\underline{q}_{y}\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy}\leq\bar{q}_{y},\
\forall y\in\mathcal{Y},$ $\displaystyle\pi_{xy}\geq 0,\
\forall\\{x,y\\}\in\mathcal{E},$

where $d_{xy}:\mathbb{R}_{+}\rightarrow\mathbb{R}$ and
$s_{xy}:\mathbb{R}_{+}\rightarrow\mathbb{R}$ are utility functions for target
node $x$ and source node $y$, respectively. They capture chosen utilities at
the UAV and target side. The upper and lower bounds are shown by
$\bar{p}_{x}$, $\bar{q}_{y}$ and $\underline{p}_{x}$, $\underline{q}_{y}$,
respectively. These bounds mean that target $x$ will not take in more than
$\bar{p}_{x}$ and take in less than $\underline{p}_{x}$. It follows similarly
for the source nodes.

###  II-B Distributed Discrete Optimal Transport

To form a distributed algorithm such that each UAV can solve its own problem,
we first make adjustments to the variables. First, the ancillary variables
$\pi_{xy,d}$ and $\pi_{xy,s}$ are introduced, where subscripts $d$ and $s$
indicate the relation to destination or source nodes, respectively. We then
set $\pi_{xy}=\pi_{xy,d}$ and $\pi_{xy,s}=\pi_{xy}$, doing so indicates that
the proposed solutions by the target and source nodes are consistent.Then, the
opposite of the original functions is taken to make a convex optimization
problem under the following assumption:

######  Assumption 1.

The utility functions $d_{xy}$ and $s_{xy}$ must be concave and monotonically
increasing.

With this assumption and the introduction of ancillary variables, we can
reformulate (1) to the following:

$\displaystyle\min_{\Pi_{t}\in\mathcal{F}_{t},\Pi_{s}\in\mathcal{F}_{s},\Pi}$
$\displaystyle-\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}d_{xy}(\pi_{xy,d})-\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}s_{xy}(\pi_{xy,s})$
(2) $\displaystyle\mathrm{s.t.}$ $\displaystyle\pi_{xy,d}=\pi_{xy},\
\forall\\{x,y\\}\in\mathcal{E},$ $\displaystyle\pi_{xy,s}=\pi_{xy},\
\forall\\{x,y\\}\in\mathcal{E},$

where $\Pi_{t}:=\\{\pi_{xy,d}\\}_{x\in\mathcal{X}_{y},y\in\mathcal{Y}}$,
$\Pi_{s}:=\\{\pi_{xy,s}\\}_{x\in\mathcal{X},y\in\mathcal{Y}_{x}}$,
$\mathcal{F}_{t}:=\\{\Pi_{t}|\pi_{xy,d}\geq
0,\underline{p}_{x}\leq\sum_{y\in\mathcal{Y}_{x}}\pi_{xy,d}\leq\bar{p}_{x},\
\\{x,y\\}\in\mathcal{E}\\}$, and $\mathcal{F}_{s}:=\\{\Pi_{s}|\pi_{xy,s}\geq
0,\underline{q}_{y}\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy,s}\leq\bar{q}_{y},\
\\{x,y\\}\in\mathcal{E}\\}$.

With the convex optimization problem in (2) we can apply the Alternating
Direction Method of Multipliers (ADMM) to form a distributed algorithm [17].

We first form a Lagrangian with $\alpha_{xy,d}$ and $\alpha_{xy,s}$ as
Lagrangian multipliers for the constraints $\pi_{xy,d}=\pi_{xy}$ and
$\pi_{xy}=\pi_{xy,s}$.

$\begin{split}&L\left(\Pi_{t},\Pi_{s},\Pi,\alpha_{xy,d},\alpha_{xy,s}\right)=-\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}d_{xy}(\pi_{xy,d})\\\
&-\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}s_{xy}(\pi_{xy,s})+\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}\alpha_{xy,d}(\pi_{xy,d}-\pi_{xy})\\\
&+\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}\alpha_{xy,s}(\pi_{xy}-\pi_{xy,s})+\frac{\eta}{2}\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}(\pi_{xy,d}-\pi_{xy})^{2}\\\
&+\frac{\eta}{2}\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}(\pi_{xy}-\pi_{xy,s})^{2},\end{split}$
(3)

where $\eta>0$ is a positive scalar controlling the convergence rate of the
following algorithm. Note that in (3), the last two terms
$\frac{\eta}{2}\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}_{x}}(\pi_{xy,d}-\pi_{xy})^{2}$
and
$\frac{\eta}{2}\sum_{y\in\mathcal{Y}}\sum_{x\in\mathcal{X}_{y}}(\pi_{xy}-\pi_{xy,s})^{2}$,
acting as penalization, are quadratic. Hence, the Lagrangian function $L$ is
strictly convex, ensuring the existence of a unique optimal solution.

Next, we apply ADMM to the minimization problem in (2) with the Lagrangian to
form the distributed algorithm:

######  Proposition 1.

The simplified iterative steps of applying ADMM to problem (2) are summarized
as follows, the simplification comes from [7]:

$\begin{split}\Pi_{x,d}(k+1)&\in\arg\min_{\Pi_{x,d}\in\mathcal{F}_{x,d}}-\sum_{y\in\mathcal{Y}_{x}}d_{xy}(\pi_{xy,d})\\\
&+\sum_{y\in\mathcal{Y}_{x}}\alpha_{xy}(k)\pi_{xy,d}+\frac{\eta}{2}\sum_{y\in\mathcal{Y}_{x}}\left(\pi_{xy,d}-\pi_{xy}(k)\right)^{2},\end{split}$
(4)
$\begin{split}\Pi_{y,s}(k+&1)\in\arg\min_{\Pi_{y,s}\in\mathcal{F}_{y,s}}-\sum_{x\in\mathcal{X}_{y}}s_{xy}(\pi_{xy,s})\\\
&-\sum_{x\in\mathcal{X}_{y}}\alpha_{xy}(k)\pi_{xy,s}+\frac{\eta}{2}\sum_{x\in\mathcal{X}_{y}}\left(\pi_{xy}(k)-\pi_{xy,s}\right)^{2},\end{split}$
(5)
$\begin{split}\pi_{xy}(k+1)=\frac{1}{2}\left(\pi_{xy,d}(k+1)+\pi_{xy,s}(k+1)\right),\end{split}$
(6)
$\begin{split}\alpha_{xy}(k+1)=\alpha_{xy}(k)+\frac{\eta}{2}\left(\pi_{xy,d}(k+1)-\pi_{xy,s}(k+1)\right),\end{split}$
(7)

where $\Pi_{\tilde{x},t}:=\\{\pi_{xy,d}\\}_{y\in\mathcal{Y}_{x},x=\tilde{x}}$
represents the solution at target node $\tilde{x}\in\mathcal{X}$, and
$\Pi_{\tilde{y},s}:=\\{\pi_{xy,s}\\}_{x\in\mathcal{X}_{y},y=\tilde{y}}$
represents the proposed solution at source node $\tilde{y}\in\mathcal{Y}$. In
addition, $\mathcal{F}_{x,d}:=\\{\Pi_{x,d}|\pi_{xy,d}\geq
0,y\in\mathcal{Y}_{x},\underline{p}_{x}\leq\sum_{y\in\mathcal{Y}_{x}}\pi_{xy,d}\leq\bar{p}_{x}\\}$,
and $\mathcal{F}_{y,s}:=\\{\Pi_{y,s}|\pi_{xy,s}\geq
0,x\in\mathcal{X}_{y},\underline{q}_{y}\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy,s}\leq\bar{q}_{y}\\}$.

For convenience, we provide these steps in the form of an algorithm:

Algorithm 1 Distributed OT Algorithm

1:while $\Pi_{x,d}$ and $\Pi_{y,s}$ not converging do

2: Compute $\Pi_{x,d}(k+1)$ using (4), for all $x\in\mathcal{X}_{y}$

3: Compute $\Pi_{y,s}(k+1)$ using (5), for all $y\in\mathcal{Y}_{x}$

4: Compute $\pi_{xy}(k+1)$ using (6), for all $\\{x,y\\}\in\mathcal{E}$

5: Compute $\alpha_{xy}(k+1)$ using (7), for all $\\{x,y\\}\in\mathcal{E}$

6:end while

7:return $\pi_{xy}(k+1)$, for all $\\{x,y\\}\in\mathcal{E}$

##  III Dynamic Distributed Optimal Transport

In this section, we adapt the above distributed OT framework for applications
to the UAV waypoint matching problem. For this application, consider the UAVs
as the source and the set of waypoints as targets where each waypoint is a
target node.

###  III-A Parameter Adjustments

To fit Algorithm 1 to the matching problem the constraints must be defined
more rigidly. First, consider that each waypoint should only have one UAV at
it at any given time, thus the constraint
$\underline{q}_{y}\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy,s}\leq\bar{q}_{y}$
becomes $0\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy,s}\leq 1$. On the source side,
every UAV must go to a waypoint thus
$\underline{p}_{x}\leq\sum_{y\in\mathcal{Y}_{x}}\pi_{xy,d}\leq\bar{p}_{x}$
becomes $\sum_{y\in\mathcal{Y}_{x}}=1$.

We also consider the network under linear parameters, making
$d_{xy}(\pi_{xy,d})=\gamma_{xy}\pi_{xy,d}$ and
$s_{xy}(\pi_{xy,s})=\delta_{xy}\pi_{xy,s}$. These parameters satisfy
Assumption 1 and are more applicable to the parameters we see in a UAV to
waypoint matching problem.

###  III-B Dynamic Optimal Transport

The network structure and parameters need to be updated often; for example,
every time a UAV reaches a the waypoint, the network parameters need to be
updated so that waypoint is not visited by any agents in the swarm again.
Also, the parameters of the nodes need to be updated based on the information
gained at the waypoint. Other events like a UAV needing to land for a battery
replacement can also affect the network. We develop a dynamic algorithm by
adding in a time step to account for this. With the dynamic algorithm, the UAV
can constantly be iterating, and when parameters are updated, it will converge
to a new optimal solution without having to start a new set of iterations.
Once the algorithm has converged, the matching scheme can be captured, and the
UAV will have its next waypoint. We introduce a new set
$t\in\mathcal{T}=\\{1,2,...,T\\}$ that indicates a time $t$. From the
modifications in subsections III-A and III-B, we obtain the following
iterative steps:

$\begin{split}\Pi_{x,d}(k+1)&\in\arg\min_{\Pi_{x,d}\in\mathcal{F}_{x,d}}-\sum_{t\in\mathcal{T}}\sum_{y\in\mathcal{Y}_{x}}\gamma_{xy}\pi_{xy,d}^{t}+\sum_{t\in\mathcal{T}}\sum_{y\in\mathcal{Y}_{x}}\\\
&\alpha_{xy}(k)\pi_{xy,d}^{t}+\frac{\eta}{2}\sum_{t\in\mathcal{T}}\sum_{y\in\mathcal{Y}_{x}}\left(\pi_{xy,d}^{t}-\pi_{xy}^{t}(k)\right)^{2},\end{split}$
(8)
$\begin{split}\Pi_{y,s}(k+&1)\in\arg\min_{\Pi_{y,s}\in\mathcal{F}_{y,s}}-\sum_{t\in\mathcal{T}}\sum_{x\in\mathcal{X}_{y}}\delta_{xy}\pi_{xy,s}^{t}-\sum_{t\in\mathcal{T}}\sum_{x\in\mathcal{X}_{y}}\\\
&\alpha_{xy}(k)\pi_{xy,s}^{t}+\frac{\eta}{2}\sum_{t\in\mathcal{T}}\sum_{x\in\mathcal{X}_{y}}\left(\pi_{xy}^{t}(k)-\pi_{xy,s}^{t}\right)^{2}\end{split}$
(9)
$\begin{split}\pi_{xy}^{t}(k+1)=\frac{1}{2}\left(\pi_{xy,d}^{t}(k+1)+\pi_{xy,s}^{t}(k+1)\right),\end{split}$
(10)
$\begin{split}\alpha_{xy}(k+1)=\alpha_{xy}(k)+\frac{\eta}{2}\left(\pi_{xy,d}^{t}(k+1)-\pi_{xy,s}^{t}(k+1)\right),\end{split}$
(11)

where
$\Pi_{\tilde{x},t}:=\\{\pi_{xy,d}\\}_{y\in\mathcal{Y}_{x},x=\tilde{x},t\in\mathcal{T}}$
represents the solution at target node $\tilde{x}\in\mathcal{X}$, and
$\Pi_{\tilde{y},s}:=\\{\pi_{xy,s}\\}_{x\in\mathcal{X}_{y},y=\tilde{y},t\in\mathcal{T}}$
represents the proposed solution at source node $\tilde{y}\in\mathcal{Y}$. In
addition, $\mathcal{F}_{x,d}:=\\{\Pi_{x,d}|\pi_{xy,d}\geq
0,y\in\mathcal{Y}_{x},\sum_{y\in\mathcal{Y}_{x}}\pi_{xy,d}^{t}=1\\}$, and
$\mathcal{F}_{y,s}:=\\{\Pi_{y,s}|\pi_{xy,s}\geq
0,x\in\mathcal{X}_{y},0\leq\sum_{x\in\mathcal{X}_{y}}\pi_{xy,s}\leq 1\\}$

Figure 2: The distributed algorithm converges to the same solution of the
centralized algorithm by showing the aggregation of the utility parameters in
the network.

Once the iterative steps have converged, the UAV will communicate what
waypoint it is heading to. When it has reached its calculated waypoint, it can
take a reading from sensors or preform a required task and communicate the
information from that point to the other UAVs. The other agents can then
update the parameters $\delta_{xy}$ connecting to the waypoint node. Once the
update is made, the iterative steps of the algorithm will begin to converge
again. We also note that the iterative steps of the algorithm will not stop
iterating. This ensures that a new optimal solution is found even if something
out of the ordinary occurs. The compromised UAV can communicate this and the
other UAVs can update their parameters and converge to a new solution rather
than having to wait until they have arrived at their waypoint to find out that
another UAV has dropped out. Also, the agents will most likely reach their
waypoints asynchronously. With a dynamic algorithm, the UAVs will not need to
wait for the other agents to reach their destinations to calculate their next
waypoint. Instead, they can calculate the optimal waypoint for themselves at
that time with the given information. For convenience, this is summarized in
Algorithm 2.

Algorithm 2 Dynamic Distributed OT Algorithm for UAVs

1:while UAV is in the air do

2: Compute $\Pi_{x,d}(k+1)$ using (8), for all $x\in\mathcal{X}_{y}$

3: Compute $\Pi_{y,s}(k+1)$ using (9), for all $y\in\mathcal{Y}_{x}$

4: Compute $\pi_{xy}(k+1)$ using (10), for all $\\{x,y\\}\in\mathcal{E}$

5: Compute $\alpha_{xy}(k+1)$ using (11), for all $\\{x,y\\}\in\mathcal{E}$

6: if Convergence is reached then

7: Interpret waypoint from $\pi_{xy}(k)$

8: Go to waypoint

9: Communicate next waypoint to swarm

10: end if

11: if At waypoint then

12: Update $\gamma_{xy}$ and $\delta_{xy}$

13: Update network structure

14: Communicate information at waypoint

15: end if

16:end while

##  IV Preliminary Case Study

In this section, we corroborate our algorithm by showing its convergence and
how it considers the parameters of the UAVs and waypoints. For example,
consider a swarm with three UAVs and ten waypoints and randomly generated
linear parameters between 1 and 10 for both $\gamma_{xy}$ and $\delta_{xy}$:

$\displaystyle\delta_{xy}=\begin{bmatrix}5&4&2&6&3&7&2&10&9&1\\\
8&2&4&5&9&5&2&4&9&2\\\ 1&1&4&7&1&6&9&7&1&9\end{bmatrix}$
$\displaystyle\gamma_{xy}=\begin{bmatrix}5&9&5&2&4&9&2&5&7&9\\\
7&1&6&9&7&1&9&10&4&1\\\ 3&7&2&10&9&1&1&6&7&8\end{bmatrix}$

We first show that the distributed algorithm converges to the solution of a
centralized algorithm. The convergence factor is set relatively high with
$\eta=10$, this allows the algorithm to converge more quickly, and time is
saved by shortening the number of iterations. The convergence of the
distributed algorithm in Alg. 1 to the centralized algorithm in (2) is shown
in Fig 2.

Once the algorithm has converged it will output a 3-by-10 matrix of zeros and
ones. The indices of the ones in the matrix indicate which UAV should go to
which waypoint. In this case, UAV 1 will be sent to waypoint 8, UAV 2 will be
sent to waypoint 5 and UAV 3 will be sent to waypoint 10.

We can corroborate the output of the algorithm by examining the aggregated
utility of the parameters $d_{xy}$ and $s_{xy}$.

$\delta_{xy}+\gamma_{xy}=\begin{bmatrix}6&10&9&14&6&12&5&17&14&3\\\
13&9&13&15&17&15&4&7&10&8\\\ 11&5&5&15&3&9&10&10&7&16\end{bmatrix}$

As this example is relatively simple, we can easily see that the algorithm is
sending the UAV to the waypoint with the highest utility as we expect it do.

###  IV-A Dynamic Algorithm

Next we verify the results of the online algorithm. In this study, at each
iteration the network is updated to exclude the waypoints that have already
been visited, this ensures that the waypoint will not be visited twice by any
of the UAVs. The parameters $\gamma_{xy}$ and $\delta_{xy}$ are also updated
to capture the new updated distances and importance of the surrounding
waypoints. For the parameter updates we randomly generate new parameters for
this study. These updates to the network are made after every 250 iterations,
but the updates to the parameters or network structure can be made at any
iteration. The algorithm converges to the new optimal solution when the
parameters are updated, as shown in Fig. 3

Figure 3: The online algorithm converges to the new optimal solution even when
parameters are adjusted while iterating.

The online algorithm is advantageous for this application because the network
will update often. Every time a UAV reaches a waypoint it will communicate
sensor readings and location to the other agents which requires the network to
update so that other UAVs can calculate their optimal waypoint. As shown in
Fig. 3 the algorithm converges in about 50 iterations, so updates can be made
often and a new optimal solution can be computed with relative ease.

##  V Simulation

(a) Waypoint Map

(b) Greedy Allocation

(c) Optimal Transport Algorithm

Figure 4: (a) shows the location of the twelve waypoints, blue indicates a
waypoint where there is no chemical present, orange indicates a waypoint with
chemical, (b) shows the path of the UAVs with the greedy algorithm, (c) shows
the path of the UAVs using Alg.2

In this section we compare the algorithm outline in Alg. 2 to a greedy
assignment algorithm in simulation. We simulate a swarm of UAVs with a
chemical sensor that will take a reading at multiple waypoints to find where
the chemical is strongest. We assume that the UAVs are altitude deconflicted
in this exercise.

###  V-A Setup

Consider an environment with three UAVs and twelve waypoints to be scanned,
three of those waypoints are considered to have a chemical. The actual
chemical sensor records a numerical reading, in simulation we report a binary
reading, meaning there is chemical present or there is not, for simplicity.
The waypoints correspond to an area with a field and a few small buildings and
the area around the buildings is the area where the chemical is considered
present, as shown in Fig. 4(a).

We compare our algorithm to a greedy assignment algorithm. The greedy
algorithm assigns four waypoints to each UAVs before they start. It assigns
waypoints solely based on distance. The point that is closest and not taken by
another UAV already will be assigned to that agent. While this algorithm does
consider distance it does not consider it at run time so a UAV maybe a
assigned the closest waypoint without the consideration of the other UAVs. The
algorithm is greedy because the agent that comes online first will choose the
waypoints that are best for itself without considering the other UAVs. Thus a
waypoint may be assigned to one UAV even if it is far closer to another UAV.

For the optimal transport algorithm in Alg. 2, the UAVs need to be assigned
their first waypoint since they all take off at the same location. The first
waypoint is assigned such that UAV one goes to waypoint one, UAV two goes to
waypoint two, and so on. We do this since the distances to the waypoints will
be the same and there is nothing known about the chemical sensor readings at
the waypoints.

###  V-B Results

The allocation schemes of the greedy and dynamic distributed optimal transport
algorithms are show in Fig. 4(b) and Fig. 4(c) respectively. Table I shows the
distances each UAV traveled between waypoints in meters.

UAV | Greedy | OT
---|---|---
1 | 255m | 108m
2 | 122m | 59m
3 | 49m | 179m
Total | 426m | 346m
TABLE I: Distances traveled

The greedy algorithm has a total distance traveled of 426 meters while the
optimal transport algorithm has a total distance traveled of 346 meters. The
smaller travel distance is essential for preserving battery life. We also
highlight that the distance are more varied even among the greedy algorithm
than with the optimal transport algorithm. The greedy algorithm has UAV one
traveling over double the distance of the other UAVs using more of that
batteries life to scan the same number of waypoints as the other agents. The
lesser total distance of the OT algorithm means that the scan of the waypoints
will be completed quicker, saving time. We note that the algorithm does take
time to converge in practice, in simulation it take around two seconds to
converge but realize this could be longer on the constrained hardware of a
UAV. These metrics show the effectiveness of the algorithm and provides a
better way to allocate the waypoints to each agent in a swarm.

##  VI Conclusion & Future Work

In this paper we developed a dynamic and distributed algorithm for the
efficient matching of UAV swarm agents to waypoints. The algorithm is capable
of considering a parameter on both the agent and waypoint side and can
converge to a new optimal matching scheme when parameters in the network
change such as a UAV needing to land for a new battery or updating parameters
based on sensor readings. The simulation shows the effectiveness of this
algorithm and showed that the developed algorithm was able to find a more
efficient way of navigating the waypoints than a greedy assignment algorithm
by cutting down the total distance that the UAVs needed to travel thus saving
battery life and shortening mission time.

Future work includes the developing a similar algorithm for a more complex use
case such as a heterogeneous swarm of robots. We also plan to cut down on time
to reach convergence making the algorithm better suited for actual hardware.
We also consider extending the security measures of the swarm especially
looking into communication between the swarm with a MPU 5 or equivalent,
extending or combining simulation with other swarm missions and building
effective yet flexible templates for mission commanders to support their
specific parameters.

##  VII Acknowledgements

This research was sponsored by the Army Research Laboratory and was
accomplished under Cooperative Agreement Number W911NF-21-2-0281. The views
and conclusions contained in this document are those of the authors and should
not be interpreted as representing the official policies, either expressed or
implied, of the Army Research Laboratory or the U.S. Government. The U.S.
Government is authorized to reproduce and distribute reprints for Government
purposes notwithstanding any copyright notation herein.

## References

  * [1] A. Sanders, “Drone swarms,” _Masters Thessis_ , 2017.
  * [2] A. Galichon, _Optimal Transport Methods in Economics_. Princeton University Press, 2018.
  * [3] R. Arnold, K. Carey, B. Abruzzo, and C. Korpela, “What is a robot swarm: a definition for swarming robotics,” in _2019 IEEE 10th Annual Ubiquitous Computing, Electronics & Mobile Communication Conference (UEMCON)_. IEEE, 2019, pp. 0074–0081.
  * [4] R. Arnold, E. Mezzacappa, M. Jablonski, J. Jablonski, and B. Abruzzo, “Performance comparison of decentralized undirected swarms versus centralized directed swarms at different levels of quality of knowledge,” in _2021 IEEE International Symposium on Technologies for Homeland Security (HST)_. IEEE, 2021, pp. 1–9.
  * [5] T. Brick, M. Lanham, A. Kopeikin, C. Korpela, and R. Morales, “Zero to swarm: Integrating suas swarming into a multi-disciplinary engineering program,” in _2018 International Conference on Unmanned Aircraft Systems (ICUAS)_ , 2018, pp. 308–314.
  * [6] E. Şahin, “Swarm robotics: From sources of inspiration to domains of application,” in _Swarm Robotics_ , E. Şahin and W. M. Spears, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2005, pp. 10–20.
  * [7] R. Zhang and Q. Zhu, “Consensus-based distributed discrete optimal transport for decentralized resource matching,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 5, no. 3, pp. 511–524, 2019.
  * [8] J. Hughes and J. Chen, “Fair and distributed dynamic optimal transport for resource allocation over networks,” in _55th Annual Conference on Information Sciences and Systems (CISS)_ , 2021.
  * [9] B. Açıkmeşe and D. S. Bayard, “Probabilistic swarm guidance for collaborative autonomous agents,” in _2014 American Control Conference_ , 2014, pp. 477–482.
  * [10] V. Krishnan and S. Martínez, “Distributed optimal transport for the deployment of swarms,” in _2018 IEEE Conference on Decision and Control (CDC)_ , 2018, pp. 4583–4588.
  * [11] Q. Wang and X. Mao, “Dynamic task allocation method of swarm robots based on optimal mass transport theory,” _Symmetry_ , vol. 12, no. 10, 2020. [Online]. Available: https://www.mdpi.com/2073-8994/12/10/1682
  * [12] B. Savidge, A. Kopeikin, R. Arnold, and D. Larkin, “Uas swarm shares survey data to expedite coordinated mapping of radiation hotspots,” in _2019 IEEE International Symposium on Technologies for Homeland Security (HST)_ , 2019, pp. 1–7.
  * [13] A. Kopeikin, S. Heider, D. Larkin, C. Korpela, R. Morales, and J. E. Bluman, _Unmanned Aircraft System Swarm for Radiological and Imagery Data Collection_. [Online]. Available: https://arc.aiaa.org/doi/abs/10.2514/6.2019-2286
  * [14] M. Ghorbanzadeh, A. Abdelhadi, and C. Clancy, _Distributed Resource Allocation_. Cham: Springer International Publishing, 2017, pp. 61–91. [Online]. Available: https://doi.org/10.1007/978-3-319-46267-7_4
  * [15] D. Niu and B. Li, “An efficient distributed algorithm for resource allocation in large-scale coupled systems,” in _2013 Proceedings IEEE INFOCOM_ , 2013, pp. 1501–1509.
  * [16] D. A. Schmidt, C. Shi, R. A. Berry, M. L. Honig, and W. Utschick, “Distributed resource allocation schemes,” _IEEE Signal Processing Magazine_ , vol. 26, no. 5, pp. 53–63, 2009.
  * [17] S. Boyd, N. Parikh, and E. Chu, _Distributed optimization and statistical learning via the alternating direction method of multipliers_. Now Publishers Inc, 2011.


* 54 Mann, A. W. _et al._ How to Constrain Your M Dwarf. II. The Mass-Luminosity-Metallicity Relation from 0.075 to 0.70 Solar Masses.  _The Astrophysical Journal_ 871, 63 (2019).  URL http://adsabs.harvard.edu/abs/2019ApJ...871...63M.
  * 55 Berger, T. A. _et al._ The Gaia-Kepler Stellar Properties Catalog. I. Homogeneous Fundamental Properties for 186,301 Kepler Stars.  _arXiv:2001.07737 [astro-ph]_ (2020).  URL http://arxiv.org/abs/2001.07737.  ArXiv: 2001.07737.
  * 56 Goodman, J. & Weare, J.  Ensemble samplers with affine invariance.  _Communications in Applied Mathematics and Computational Science_ 5, 65–80 (2010).  URL https://msp.org/camcos/2010/5-1/p04.xhtml.  Publisher: Mathematical Sciences Publishers.
  * 57 Nelson, B. E., Ford, E. B. & Payne, M. J.  RUN DMC: An efficient, parallel code for analyzing Radial Velocity Observations using N-body Integrations and Differential Evolution Markov chain Monte Carlo.  _The Astrophysical Journal Supplement Series_ 210, 11 (2013).  URL http://arxiv.org/abs/1311.5229.  ArXiv:1311.5229 [astro-ph].
  * 58 Heyl, J. S. & Gladman, B. J.  Using long-term transit timing to detect terrestrial planets.  _Monthly Notices of the Royal Astronomical Society_ 377, 1511–1519 (2007).  URL http://adsabs.harvard.edu/abs/2007MNRAS.377.1511H.
  * 59 Rein, H. & Liu, S.-F.  REBOUND: an open-source multi-purpose N-body code for collisional dynamics.  _Astronomy and Astrophysics_ 537, A128 (2012).  URL http://adsabs.harvard.edu/abs/2012A%26A...537A.128R.
  * 60 Rein, H. & Tamayo, D.  WHFAST: a fast and unbiased implementation of a symplectic Wisdom-Holman integrator for long-term gravitational simulations.  _Monthly Notices of the Royal Astronomical Society_ 452, 376–388 (2015).  URL http://adsabs.harvard.edu/abs/2015MNRAS.452..376R.
  * 61 Mandel, K. & Agol, E.  Analytic Light Curves for Planetary Transit Searches.  _The Astrophysical Journal Letters_ 580, L171–L175 (2002).  URL http://adsabs.harvard.edu/abs/2002ApJ...580L.171M.
  * 62 Wang, J., Fischer, D. A., Xie, J.-W. & Ciardi, D. R.  Influence of Stellar Multiplicity on Planet Formation. IV. Adaptive Optics Imaging of Kepler Stars with Multiple Transiting Planet Candidates.  _The Astrophysical Journal_ 813, 130 (2015).  URL https://ui.adsabs.harvard.edu/abs/2015ApJ...813..130W.  ADS Bibcode: 2015ApJ…813..130W.
  * 63 Vogt, S. S. _et al._ HIRES: the high-resolution echelle spectrometer on the Keck 10-m Telescope 2198, 362 (1994).  URL http://adsabs.harvard.edu/abs/1994SPIE.2198..362V.  Conference Name: Instrumentation in Astronomy VIII.
  * 64 Butler, R. P. _et al._ Attaining Doppler Precision of 3 M s-1.  _Publications of the Astronomical Society of the Pacific_ 108, 500 (1996).  URL http://adsabs.harvard.edu/abs/1996PASP..108..500B.
  * 65 Fulton, B. J., Petigura, E. A., Blunt, S. & Sinukoff, E.  RadVel: The Radial Velocity Modeling Toolkit.  _Publications of the Astronomical Society of the Pacific_ 130, 044504 (2018).  URL https://ui.adsabs.harvard.edu/abs/2018PASP..130d4504F/abstract.
  * 66 Amado, P. J. _et al._ The CARMENES search for exoplanets around M dwarfs. Two terrestrial planets orbiting G 264-012 and one terrestrial planet orbiting Gl 393.  _arXiv e-prints_ 2105, arXiv:2105.13785 (2021).  URL http://adsabs.harvard.edu/abs/2021arXiv210513785A.
  * 67 Ahrer, E. _et al._ The HARPS search for southern extra-solar planets XLV. Two Neptune mass planets orbiting HD 13808: a study of stellar activity modelling’s impact on planet detection.  _arXiv e-prints_ 2102, arXiv:2102.03387 (2021).  URL http://adsabs.harvard.edu/abs/2021arXiv210203387A.
  * 68 Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W. & O’Neil, M.  Fast Direct Methods for Gaussian Processes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 38 (2015).  URL http://adsabs.harvard.edu/abs/2015ITPAM..38..252A.
  * 69 McQuillan, A., Aigrain, S. & Mazeh, T.  Measuring the rotation period distribution of field M dwarfs with Kepler.  _Monthly Notices of the Royal Astronomical Society_ 432, 1203–1216 (2013).  URL http://adsabs.harvard.edu/abs/2013MNRAS.432.1203M.
  * 70 McQuillan, A., Mazeh, T. & Aigrain, S.  Stellar Rotation Periods of the Kepler Objects of Interest: A Dearth of Close-in Planets around Fast Rotators.  _The Astrophysical Journal Letters_ 775, L11 (2013).  URL http://adsabs.harvard.edu/abs/2013ApJ...775L..11M.
  * 71 Benneke, B.  Strict Upper Limits on the Carbon-to-Oxygen Ratios of Eight Hot Jupiters from Self-Consistent Atmospheric Retrieval.  _arXiv e-prints_ 1504, arXiv:1504.07655 (2015).  URL http://adsabs.harvard.edu/abs/2015arXiv150407655B.
  * 72 Benneke, B. & Seager, S.  Atmospheric Retrieval for Super-Earths: Uniquely Constraining the Atmospheric Composition with Transmission Spectroscopy.  _The Astrophysical Journal_ 753, 100 (2012).  URL http://adsabs.harvard.edu/abs/2012ApJ...753..100B.
  * 73 Benneke, B. & Seager, S.  How to Distinguish between Cloudy Mini-Neptunes and Water/Volatile-dominated Super-Earths.  _The Astrophysical Journal_ 778, 153 (2013).  URL http://adsabs.harvard.edu/abs/2013ApJ...778..153B.
  * 74 Feroz, F., Hobson, M. P. & Bridges, M.  MULTINEST: an efficient and robust Bayesian inference tool for cosmology and particle physics.  _Monthly Notices of the Royal Astronomical Society_ 398, 1601–1614 (2009).  URL http://adsabs.harvard.edu/abs/2009MNRAS.398.1601F.
  * 75 Shaw, J. R., Bridges, M. & Hobson, M. P.  Efficient Bayesian inference for multimodal problems in cosmology.  _Monthly Notices of the Royal Astronomical Society_ 378, 1365–1370 (2007).  URL https://ui.adsabs.harvard.edu/abs/2007MNRAS.378.1365S.
  * 76 Mukherjee, P., Parkinson, D. & Liddle, A. R.  A Nested Sampling Algorithm for Cosmological Model Selection.  _The Astrophysical Journal_ 638, L51–L54 (2006).  URL https://ui.adsabs.harvard.edu/abs/2006ApJ...638L..51M.
  * 77 Skilling, J.  Nested Sampling 735, 395–405 (2004).  URL http://adsabs.harvard.edu/abs/2004AIPC..735..395S.  Conference Name: Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering.
  * 78 Trotta, R.  Bayes in the sky: Bayesian inference and model selection in cosmology.  _Contemporary Physics_ 49, 71 (2008).  URL https://ui.adsabs.harvard.edu/abs/2008ConPh..49...71T/abstract.
  * 79 Line, M. R. & Parmentier, V.  The Influence of Nonuniform Cloud Cover on Transit Transmission Spectra.  _The Astrophysical Journal_ 820, 78 (2016).  URL http://adsabs.harvard.edu/abs/2016ApJ...820...78L.
  * 80 Miller-Ricci, E., Seager, S. & Sasselov, D.  The Atmospheric Signatures of Super-Earths: How to Distinguish Between Hydrogen-Rich and Hydrogen-Poor Atmospheres.  _The Astrophysical Journal_ 690, 1056–1067 (2009).  URL http://adsabs.harvard.edu/abs/2009ApJ...690.1056M.
  * 81 Thorngren, D., Gao, P. & Fortney, J. J.  The Intrinsic Temperature and Radiative-Convective Boundary Depth in the Atmospheres of Hot Jupiters.  _The Astrophysical Journal Letters_ 884, L6 (2019).  URL http://adsabs.harvard.edu/abs/2019ApJ...884L...6T.
  * 82 Chabrier, G., Mazevet, S. & Soubiran, F.  A New Equation of State for Dense Hydrogen-Helium Mixtures.  _The Astrophysical Journal_ 872, 51 (2019).  URL http://adsabs.harvard.edu/abs/2019ApJ...872...51C.
  * 83 Mazevet, S., Licari, A., Chabrier, G. & Potekhin, A. Y.  Ab initio based equation of state of dense water for planetary and exoplanetary modeling.  _Astronomy and Astrophysics_ 621, A128 (2019).  URL http://adsabs.harvard.edu/abs/2019A%26A...621A.128M.
  * 84 Valencia, D., Guillot, T., Parmentier, V. & Freedman, R. S.  BULK COMPOSITION OF GJ 1214b AND OTHER SUB-NEPTUNE EXOPLANETS.  _The Astrophysical Journal_ 775, 10 (2013).  URL https://iopscience.iop.org/article/10.1088/0004-637X/775/1/10.
  * 85 Madhusudhan, N., Nixon, M. C., Welbanks, L., Piette, A. A. A. & Booth, R. A.  The Interior and Atmosphere of the Habitable-zone Exoplanet K2-18b.  _The Astrophysical Journal Letters_ 891, L7 (2020).  URL http://adsabs.harvard.edu/abs/2020ApJ...891L...7M.
  * 86 Lopez, E. D. & Fortney, J. J.  Understanding the Mass-Radius Relation for Sub-neptunes: Radius as a Proxy for Composition.  _The Astrophysical Journal_ 792, 1 (2014).  URL http://adsabs.harvard.edu/abs/2014ApJ...792....1L.
  * 87 Hubbard, W. B. _et al._ Theory of Extrasolar Giant Planet Transits.  _The Astrophysical Journal_ 560, 413–419 (2001).  URL https://ui.adsabs.harvard.edu/abs/2001ApJ...560..413H.  ADS Bibcode: 2001ApJ…560..413H.
  * 88 Otegi, J. F. _et al._ Impact of the measured parameters of exoplanets on the inferred internal structure.  _Astronomy & Astrophysics, Volume 640, id.A135, <NUMPAGES>12</NUMPAGES> pp._ 640, A135 (2020).  URL https://ui.adsabs.harvard.edu/abs/2020A%26A...640A.135O/abstract.
  * 89 Lozovsky, M., Helled, R., Dorn, C. & Venturini, J.  Threshold Radii of Volatile-rich Planets.  _The Astrophysical Journal_ 866, 49 (2018).  URL https://iopscience.iop.org/article/10.3847/1538-4357/aadd09.
  * 90 Turbet, M., Ehrenreich, D., Lovis, C., Bolmont, E. & Fauchez, T.  The runaway greenhouse radius inflation effect. An observational diagnostic to probe water on Earth-sized planets and test the habitable zone concept.  _Astronomy and Astrophysics_ 628, A12 (2019).  URL http://adsabs.harvard.edu/abs/2019A%26A...628A..12T.
  * 91 Turbet, M. _et al._ Revised mass-radius relationships for water-rich rocky planets more irradiated than the runaway greenhouse limit.  _Astronomy & Astrophysics, Volume 638, id.A41, <NUMPAGES>10</NUMPAGES> pp._ 638, A41 (2020).  URL https://ui.adsabs.harvard.edu/abs/2020A%26A...638A..41T/abstract.
  * 92 Zeng, L., Sasselov, D. D. & Jacobsen, S. B.  Mass-Radius Relation for Rocky Planets Based on PREM.  _The Astrophysical Journal_ 819, 127 (2016).  URL http://adsabs.harvard.edu/abs/2016ApJ...819..127Z.
  * 93 Madhusudhan, N., Piette, A. A. A. & Constantinou, S.  Habitability and Biosignatures of Hycean Worlds.  _The Astrophysical Journal_ 918, 1 (2021).  URL https://ui.adsabs.harvard.edu/abs/2021ApJ...918....1M.  ADS Bibcode: 2021ApJ…918….1M.
  * 94 Kosiarek, M. R. _et al._ Physical Parameters of the Multi-Planet Systems HD 106315 and GJ 9827.  _arXiv:2009.03398 [astro-ph]_ (2020).  URL http://arxiv.org/abs/2009.03398.  ArXiv: 2009.03398.
  * 95 Curtis, J. L. _et al._ When Do Stalled Stars Resume Spinning Down? Advancing Gyrochronology with Ruprecht 147.  _The Astrophysical Journal_ 904, 140 (2020).  URL https://ui.adsabs.harvard.edu/abs/2020ApJ...904..140C.  ADS Bibcode: 2020ApJ…904..140C.
  * 96 Muirhead, P. S. _et al._ Characterizing the Cool Kepler Objects of Interests. New Effective Temperatures, Metallicities, Masses, and Radii of Low-mass Kepler Planet-candidate Host Stars.  _The Astrophysical Journal_ 750, L37 (2012).  URL https://ui.adsabs.harvard.edu/2012ApJ...750L..37M/abstract.
  * 97 Watson, A. J., Donahue, T. M. & Walker, J. C. G.  The dynamics of a rapidly escaping atmosphere: Applications to the evolution of Earth and Venus.  _Icarus_ 48, 150–166 (1981).  URL https://ui.adsabs.harvard.edu/abs/1981Icar...48..150W.  ADS Bibcode: 1981Icar…48..150W.
  * 98 Owen, J. E. & Wu, Y.  The Evaporation Valley in the Kepler Planets.  _The Astrophysical Journal_ 847, 29 (2017).  URL http://adsabs.harvard.edu/abs/2017ApJ...847...29O.
  * 99 Feinstein, A. D. _et al._ Flare Statistics for Young Stars from a Convolutional Neural Network Analysis of $\textit{TESS}$ Data.  _The Astronomical Journal_ 160, 219 (2020).  URL http://arxiv.org/abs/2005.07710.  ArXiv: 2005.07710.
  * 100 Ribas, I., Guinan, E. F., Güdel, M. & Audard, M.  Evolution of the Solar Activity over Time and Effects on Planetary Atmospheres. I. High-Energy Irradiances (1-1700 Å).  _The Astrophysical Journal_ 622, 680–694 (2005).  URL http://adsabs.harvard.edu/abs/2005ApJ...622..680R.
  * 101 Jackson, A. P., Davis, T. A. & Wheatley, P. J.  The coronal X-ray-age relation and its implications for the evaporation of exoplanets.  _Monthly Notices of the Royal Astronomical Society_ 422, 2024–2043 (2012).  URL http://adsabs.harvard.edu/abs/2012MNRAS.422.2024J.
  * 102 Tu, L., Johnstone, C. P., Güdel, M. & Lammer, H.  The extreme ultraviolet and X-ray Sun in Time: High-energy evolutionary tracks of a solar-like star.  _Astronomy and Astrophysics_ 577, L3 (2015).  URL http://adsabs.harvard.edu/abs/2015A%26A...577L...3T.
  * 103 Güdel, M., Guinan, E. F. & Skinner, S. L.  The X-Ray Sun in Time: A Study of the Long-Term Evolution of Coronae of Solar-Type Stars.  _The Astrophysical Journal_ 483, 947–960 (1997).  URL http://adsabs.harvard.edu/abs/1997ApJ...483..947G.
  * 104 Murray-Clay, R. A., Chiang, E. I. & Murray, N.  Atmospheric Escape From Hot Jupiters.  _The Astrophysical Journal_ 693, 23–42 (2009).  URL http://adsabs.harvard.edu/abs/2009ApJ...693...23M.
  * 105 Owen, J. E. & Jackson, A. P.  Planetary evaporation by UV & X-ray radiation: basic hydrodynamics.  _Monthly Notices of the Royal Astronomical Society_ 425, 2931–2947 (2012).  URL http://adsabs.harvard.edu/abs/2012MNRAS.425.2931O.
  * 106 Owen, J. E. & Alvarez, M. A.  UV Driven Evaporation of Close-in Planets: Energy-limited, Recombination-limited, and Photon-limited Flows.  _The Astrophysical Journal_ 816, 34 (2016).  URL http://adsabs.harvard.edu/abs/2016ApJ...816...34O.
  * 107 Erkaev, N. V. _et al._ EUV-driven mass-loss of protoplanetary cores with hydrogen-dominated atmospheres: the influences of ionization and orbital distance.  _Monthly Notices of the Royal Astronomical Society_ 460, 1300–1309 (2016).  URL https://ui.adsabs.harvard.edu/abs/2016MNRAS.460.1300E.  ADS Bibcode: 2016MNRAS.460.1300E.
  * 108 Erkaev, N. V. _et al._ Roche lobe effects on the atmospheric loss from “Hot Jupiters”.  _Astronomy and Astrophysics_ 472, 329–334 (2007).  URL http://adsabs.harvard.edu/abs/2007A%26A...472..329E.
  * 109 Johnstone, C. P., Bartel, M. & Güdel, M.  The active lives of stars: A complete description of the rotation and XUV evolution of F, G, K, and M dwarfs.  _Astronomy and Astrophysics_ 649, A96 (2021).  URL https://ui.adsabs.harvard.edu/abs/2021A&A...649A..96J/abstract.
  * 110 Schaefer, L. & Fegley, B.  Chemistry of atmospheres formed during accretion of the Earth and other terrestrial planets.  _Icarus_ 208, 438–448 (2010).  URL http://adsabs.harvard.edu/abs/2010Icar..208..438S.
  * 111 Lichtenberg, T. _et al._ Vertically resolved magma ocean-protoatmosphere evolution: H$_2$, H$_2$O, CO$_2$, CH$_4$, CO, O$_2$, and N$_2$ as primary absorbers.  _Journal of Geophysical Research: Planets_ (2021).  URL http://arxiv.org/abs/2101.10991.  ArXiv: 2101.10991.
  * 112 Sossi, P. A.  Atmospheres in the baking.  _Nature Astronomy_ 5, 535–536 (2021).  URL https://ui.adsabs.harvard.edu/abs/2021NatAs...5..535S.  ADS Bibcode: 2021NatAs…5..535S.
  * 113 Andrault, D., Monteux, J., Le Bars, M. & Samuel, H.  The deep Earth may not be cooling down.  _Earth and Planetary Science Letters, Volume 443, p. 195-203._ 443, 195 (2016).  URL https://ui.adsabs.harvard.edu/abs/2016E%26PSL.443..195A/abstract.
  * 114 Rackham, B. V., Apai, D. & Giampapa, M. S.  The Transit Light Source Effect: False Spectral Features and Incorrect Densities for M-dwarf Transiting Planets.  _The Astrophysical Journal_ 853, 122 (2018).  URL http://adsabs.harvard.edu/abs/2018ApJ...853..122R.
  * 115 Husser, T.-O. _et al._ A new extensive library of PHOENIX stellar atmospheres and synthetic spectra.  _Astronomy & Astrophysics, Volume 553, id.A6, <NUMPAGES>9</NUMPAGES> pp._ 553, A6 (2013).  URL https://ui.adsabs.harvard.edu/abs/2013A%26A...553A...6H/abstract.
  * 116 Gao, P. & Zhang, X.  Deflating Super-puffs: Impact of Photochemical Hazes on the Observed Mass-Radius Relationship of Low-mass Planets.  _The Astrophysical Journal_ 890, 93 (2020).  URL http://adsabs.harvard.edu/abs/2020ApJ...890...93G.
  * 117 Lavvas, P., Koskinen, T., Steinrueck, M. E., García Muñoz, A. & Showman, A. P.  Photochemical Hazes in Sub-Neptunian Atmospheres with a Focus on GJ 1214b.  _The Astrophysical Journal_ 878, 118 (2019).  URL http://adsabs.harvard.edu/abs/2019ApJ...878..118L.
  * 118 Piro, A. L.  Can Rocky Exoplanets with Rings Pose as Sub-Neptunes?  _The Astronomical Journal_ 156, 80 (2018).  URL http://adsabs.harvard.edu/abs/2018AJ....156...80P.
  * 119 Piro, A. L. & Vissapragada, S.  Exploring Whether Super-Puffs Can Be Explained as Ringed Exoplanets.  _The Astronomical Journal_ 159, 131 (2020).  URL http://arxiv.org/abs/1911.09673.  ArXiv: 1911.09673.
  * 120 Clausen, N. & Tilgner, A.  Dissipation in rocky planets for strong tidal forcing.  _Astronomy and Astrophysics_ 584, A60 (2015).  URL http://adsabs.harvard.edu/abs/2015A%26A...584A..60C.
  * 121 Chandrasekhar, S.  Ellipsoidal figures of equilibrium.  _The Silliman Foundation Lectures, New Haven: Yale University Press, 1969_ (1969).  URL http://adsabs.harvard.edu/abs/1969efe..book.....C.
  * 122 Tremaine, S., Touma, J. & Namouni, F.  Satellite Dynamics on the Laplace Surface.  _The Astronomical Journal_ 137, 3706–3717 (2009).  URL http://adsabs.harvard.edu/abs/2009AJ....137.3706T.
  * 123 Schlichting, H. E. & Chang, P.  Warm Saturns: On the Nature of Rings around Extrasolar Planets That Reside inside the Ice Line.  _The Astrophysical Journal_ 734, 117 (2011).  URL http://adsabs.harvard.edu/abs/2011ApJ...734..117S.
  * 124 Astropy Collaboration _et al._ Astropy: A community Python package for astronomy.  _A &A_ 558, A33 (2013).  1307.6212.
  * 125 Astropy Collaboration _et al._ The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package.  _AJ_ 156, 123 (2018).  1801.02634.
  * 126 Harris, C. R. _et al._ Array programming with NumPy.  _Nature_ 585, 357–362 (2020).  URL https://doi.org/10.1038/s41586-020-2649-2.
  * 127 Pérez, F. & Granger, B. E.  Ipython: a system for interactive scientific computing (2007).
  * 128 Hunter, J. D.  Matplotlib: A 2d graphics environment.  _Computing in Science & Engineering_ 9, 90–95 (2007).
  * 129 Tange, O.  Gnu parallel 20200722 (’privacy shield’) (2020).  URL https://doi.org/10.5281/zenodo.3956817.  GNU Parallel is a general parallelizer to run multiple serial command line programs in parallel without changing them.
  * 130 Zeng, L. & Sasselov, D.  A Detailed Model Grid for Solid Planets from 0.1 through 100 Earth Masses.  _Publications of the Astronomical Society of the Pacific_ 125, 227 (2013).  URL http://adsabs.harvard.edu/abs/2013PASP..125..227Z.
  * 131 Marcus, R. A., Sasselov, D., Hernquist, L. & Stewart, S. T.  Minimum Radii of Super-Earths: Constraints from Giant Impacts.  _The Astrophysical Journal Letters_ 712, L73–L76 (2010).  URL http://adsabs.harvard.edu/abs/2010ApJ...712L..73M.
  * 132 Lu, Y. L., Angus, R., Curtis, J. L., David, T. J. & Kiman, R.  Gyro-kinematic Ages for around 30,000 Kepler Stars.  _The Astronomical Journal_ 161, 189 (2021).  URL https://ui.adsabs.harvard.edu/abs/2021AJ....161..189L.  ADS Bibcode: 2021AJ….161..189L.
  * 133 Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G. & Andrae, R.  Estimating Distance from Parallaxes. IV. Distances to 1.33 Billion Stars in Gaia Data Release 2.  _The Astronomical Journal_ 156, 58 (2018).  URL https://ui.adsabs.harvard.edu/abs/2018AJ....156...58B.