# Photometric Determination of Main-Sequence Binaries with Gaia

A. L. Wallace${}^{\href https://orcid.org/0000-0002-6591-5290\,1}$
1School of Physics & Astronomy, Monash University, Victoria, Australia
E-mail<EMAIL_ADDRESS>
###### Abstract

Since its launch in 2013, the Gaia space telescope has provided precise
measurements of the positions and magnitudes of over 1 billion stars. This has
enabled extensive searches for stellar and sub-stellar companions through
astrometric and radial velocity measurements. However, these surveys require a
prior knowledge of any unresolved companion affecting the results which can be
identified using photometry. In this work, Gaia’s magnitude measurements are
combined with near-infrared observations from 2MASS and WISE and simulation-
based inference is applied to constrain astrophysical parameters and search
for hidden companions. This method is first tested on simulated sets of binary
stars before expanding to Gaia’s non-single star catalogue. Using this test, a
region is identified on the H-R diagram in which the method is the most
accurate and all Gaia sources within that region are analysed. This analysis
reproduces a known anti-correlation between metallicity and binary fraction.
Finally, the method is applied to the nearby star cluster M67 and, using
previous studies of the metallicity distribution, it is possible to improve
constraints on binary fraction. From this the binary fraction in the cluster
is calculated to vary from 30% in the outer cluster to 45% near the core. This
is found to be significantly higher the 23% binary fraction calculated for the
wider stellar neighbourhood.

###### keywords:

binaries: close – stars: general – stars: fundamental parameters

††pubyear: 2023††pagerange: Photometric Determination of Main-Sequence
Binaries with Gaia–Photometric Determination of Main-Sequence Binaries with
Gaia

##  1 Introduction

It has long been understood that a significant fraction of stars in our galaxy
form in multiple systems. Previous studies show about 10% of field stars are
hosting at least one companion, a fraction which more than doubles in star
clusters and star-forming regions (Duchêne et al., 2018). This substantial
difference in binary fraction between different stellar populations is of
particular interest when studying stellar evolution. By studying stellar
multiplicity in a wide range of environments such as star clusters and the
Galactic field, and taking dynamical interactions into account, this can help
test theories of star formation (Duchêne & Kraus, 2013). For instance, it has
been established that wide binaries (>100 AU) are rare in young star clusters
(Deacon & Kraus, 2020) leading to different theories of stellar interactions
and evolution to explain this discrepancy.
Studies of tight binary systems (<100 AU) are important for our understanding
both stellar evolution and planet formation (Saleh & Rasio, 2009). It is
understood that tight binary fraction, like total binary fraction, increases
with the mass of the primary star (Clark et al., 2012) and that the chemical
properties of a star may play an even larger role in determining stellar
multiplicity (Badenes et al., 2018; Mazzola et al., 2020). Due to the many
factors such as a star’s mass, composition and neighbourhood conditions
affecting the binary fraction, it is essential to form a highly accurate
picture of binary fraction as a function of these properties. Studies of tight
binaries has only been made possible in recent decades by advances in
spectroscopy and radial velocity measurements. However, with highly accurate
photometry and astrometry, it is also possible to determine stellar
multiplicity from a star’s observed colour and magnitude.
Since the launch of the Gaia Space Telescope, we now have precise measurements
of the positions of over 1 billion stars providing a detailed 3-dimensional
map of the Milky Way. The third data release (DR3), based on the first 34
months of observations (Gaia Collaboration et al., 2021), contains a wealth of
new information on existing sources including precise radial velocity (RV) and
astrometric observations. These observations are useful in searching for
stellar companions and have successfully identified previously unknown
multiple star systems (Mugrauer et al., 2022).
As well as RV and astrometry, Gaia has also provided us with highly accurate
magnitudes and colours of stars in the Milky Way. Using these magnitudes, it
is possible to identify unresolved stellar-mass companions to main-sequence
stars. Stellar multiplicity on the main-sequence can be observed on the H-R
diagram, as shown in Figure 1, for three different primary masses. The
isochrone for a single atar of solar age and metallicity is shown for
comparison.

As shown in Figure 1, as the binary mass ratio increases, the star becomes
redder and slightly brighter. Above a mass ratio of $\sim$0.7, the star’s
colour shifts to shorter wavelengths and the brightness increases rapidly,
with an equal binary having a magnitude $\sim$0.75 brighter than a single
star.
Unresolved stellar companions are an important factor to consider when
conducting astrometric surveys with Gaia as the presence of an additional
high-mass companion could hinder the search for lower mass companions if not
taken into account. By independently investigating binarity early on, we can
remove these potential obstacles from future astrometric surveys. However, a
major challenge in identifying binaries from photometry alone is the multiple
astrophysical parameters, such as metallicity, primary mass and age, all
affecting a star’s position on the H-R diagram. Previous studies have
constrained stellar properties of binaries (Traven et al., 2020) but these
were binaries that had already been identified through other methods. Other
studies have investigated how photometry can determine binarity from
photometry (Hartigan & Kenyon, 2001) but did not constrain other stellar
parameters.
In this study, simulation-based inference is applied to photometric
measurements from Gaia to constrain mass ratios of potential binaries as well
as the primary mass, metallicity and age. Degeneracies are minimised by also
including near infrared photometry from the 2MASS and WISE surveys. For the
purposes of this study, the sample is restricted to the main sequence due to
higher numbers of potential targets and more constrained stellar models for
this evolutionary stage. In the following sections, the method is introduced
and first tested on simulated data before being applied to a selection of
nearby star clusters and a large sample across the entire sky from the Gaia
DR3 catalogue to paint a clearer picture of stellar multiplicity in the local
neighbourhood.

Figure 1: Positions of binaries on the H-R diagram for primary masses of 0.6,
0.8 and 1.0 M${\odot}$. The isochrone for single main-sequence stars of solar
age (4.6 Gyr) and metallicity ([Fe/H]=0) is also shown.

##  2 Method

This section describes the method in which astrophysical parameters (primary
mass, mass ratio, metallicity, age and distance) are calculated based on
observed magnitudes. As multiple parameters are recovered from multiple
magnitudes, astrophysical parameters are denoted by the vector

$\vec{\theta}=\left[M_{1},q,\tau,\left[\mathrm{Fe/H}\right],d\right]$ (1)

where $M_{1}$ is the primary mass, $q$ is the mass ratio (always between 0 and
1), $\tau$ is the age of the system in Gyr, [Fe/H] is the metallicity relative
to the Sun and $d$ is the distance in pc. We denote the observed quantities
with the vector

$\displaystyle\vec{x}=[m_{B},m_{G},m_{R},m_{J},m_{H},m_{K},m_{W1},m_{W2},\varpi,$
(2)
$\displaystyle\sigma_{B},\sigma_{G},\sigma_{R},\sigma_{J},\sigma_{H},\sigma_{K},\sigma_{W1},\sigma_{W2},\sigma_{\varpi}]$

where the $m$’s are the apparent magnitudes in the Bp, G, Rp, J, H, Ks, W1 and
W2 wavelength bands, $\varpi$ is the parallax in mas measured by Gaia and the
$\sigma$’s are the associated errors in these observations.

In order to recover astrophysical parameters $\vec{\theta}$ from observations
$\vec{x}$, there are several approaches we could take, which rely on
possessing an accurate model of $\vec{x}$ as a function of $\vec{\theta}$. The
simplest method would be to construct a likelihood function which models the
probability a certain set of parameters produce the observed magnitudes. A set
of parameters which produces the maximum value of this probability would then
be taken as the true parameters. However, the relation between $\vec{\theta}$
and $\vec{x}$ is highly degenerate (e.g. a metal-rich single star could have
similar magnitudes to a metal-poor binary) which produces local maxima and
minima in the likelihood function. When scanning the likelihood function there
is then a risk of getting stuck on a local maximum with parameters far away
from the true $\vec{\theta}$. In order for this method to work, an initial
guess for $\vec{\theta}$ is required to be implausibly close to the true
value.

The aim of this study is to obtain the posterior distributions of the
astrophysical parameters $p(\vec{\theta}|\vec{x})$ which is a measure of the
probability of a set of parameters $\vec{\theta}$ for given data $\vec{x}$.
This distribution is generally expressed by Bayes’ Theorem:

$p(\vec{\theta}|\vec{x})=p(\vec{x}|\vec{\theta})p(\vec{theta})$ (3)

where $p(\vec{x}|\vec{\theta})$ is the probability of the data $\vec{x}$ for
given parameters $\vec{\theta}$ (the likelihood) and $p(\vec{theta})$ is the
prior or expected distribution of the parameters. This study relies on
simulation-based inference (SBI) (Cranmer et al., 2020), also known as
likelihood-free inference which aims to approximate the posterior distribution
with no assumptions about the form of the likelihood function. SBI works by
using a set of parameters $\vec{\theta}$ to generate mock data $\vec{x}$ and
training a neural network to determine a relation between $\vec{\theta}$ and
$\vec{x}$. This relation is then used to calculate an $N$-dimensional
probability function for a given set of observations, where $N$ is the number
of parameters. This method is faster than a likelihood maximisation as we rely
on an overall distribution rather than a point-by-point comparison and, as we
can estimate the density in parameter space, we can concentrate on the most
likely values for $\vec{\theta}$.

###  2.1 Simulated Training Set

Before analysing real data, a simulated set of stars is constructed with given
$M_{1}$, $q$, $\tau$, [Fe/H] and $d$. This is the set of astrophysical
parameters denoted $\vec{\theta}$ (Equation 1). For this ‘training’ set,
100,000 systems are simulated, with M1 ranging from 0.4–5 M⊙, $q$ from 0–1,
$\tau$ from 0.2–10 Gyr and [Fe/H] from -2 to +0.5.

The values for $M_{1}$ and [Fe/H] in the training set were sampled from beta
distributions designed to approximate the initial mass function and an
expected metallicity distribution which is maximised near solar metallicity.
This distribution has a probability density function given by:

$p(x;\alpha,\beta)\propto x^{\alpha-1}(1-x)^{\beta-1}$ (4)

where x is the parameter (either $M_{1}$ or [Fe/H]). The mass distribution has
$\alpha=1$, $\beta=5$ and the metallicity distribution has $\alpha=10$,
$\beta=2$. The values of $q$ and $\tau$ are sampled from uniform
distributions.

These stars are simulated with the isochrones package (Morton, 2015) to
convert the astrophysical parameters $\vec{\theta}$ into a set of absolute
magnitudes in the Gaia (Bp,G,Rp), 2MASS (J,H,Ks) and WISE (W1,W2) bands. If a
star produces an invalid result (e.g. it is too old for its mass) it is
discarded and new stars are simulated until 100,000 valid results are
produced. This trainingn set initially ignores distance and only simulates
absolute magnitudes in all bands. The training set of 100,000 stars are shown
as a histrogram on the H-R diagram in Figure 2.

Figure 2: Histogram of the simulated training set on the H-R diagram before
the addition of errors.

After producing the sample set of absolute magnitudes. Each star is assigned
100 distances between 1–1000 pc sampled from a uniform distribution. The
absolute magnitudes are then converted to apparent magnitudes. The training
now consists of 10 million stars, each with a parallax and apparent magnitudes
in eight bands.

####  2.1.1 Uncertainties

Uncertainties in parallax and magnitude are handled using a similar method to
Hahn & Melchior (2022), in which uncertainties are assigned to the simulated
data and treat these as extra variables when running the SBI. These
uncertainties are calculated using real observed uncertainties from Gaia,
2MASS and WISE. These observations give an indication of how the uncertainties
should vary as a function of magnitude. The distributions for the real data
are shown in Figure 3.

Figure 3: Distributions of magnitudes and associated errors, shown on a log
density scale from dark blue (lowest) to yellow (highest). Note parallax error
is shown as a function of G magnitude rather than parallax as the magnitude is
a key factor determining astrometric errors (and therefore errors in parallax)

(Lindegren et al., 2021).

Using this real data, a grid of colours (Bp-Rp) and apparent G magnitude is
constructed. Each simulated source is sorted into a square on this grid based
on its simulated colour and magnitude. A real source is randomly selected from
this square and its errors in parallax and magnitude are applied to the
simulated star. This ensures that the simulated source will have a set of
errors that is consistent with real data, rather than a set of independently
generate random errors.

In total, this simulation produces 18 ‘observables’ for each star (as shown in
Equation 2): a parallax, apparent magnitudes in Gaia bands ($Bp,G,Rp$), 2MASS
bands ($J,H,Ks$) and WISE bands ($W_{1},W_{2}$) and their associated
uncertainties. The set of observables is denoted $\vec{x}$.

###  2.2 Parameter Inference

The set of astrophysical parameters $\vec{\theta}$ which map to observables
$\vec{x}$ are then used to train a neural network which is used to derive the
full parameter space and approximate the posterior distribution
$p(\vec{\theta}|\vec{x})$ where $\vec{x}$ is a set of observations, either
real or simulated.
An example corner plot of the posteriors for a simulated binary is shown in
Figure 4.

Figure 4: Example corner plot showing parameter distributions for a simulated
binary star and correlation between parameters. In this example, a set of
observations are calculated using isochrones for a set of given astrophysical
parameters. The SBI method is then applied to these observations in an attempt
to reproduce the simulated parameters. The simulated ‘true’ values are then
shown for comparison.

As shown in Figure 4, the SBI method in this particular case has succeeded in
accurately recovering the correct masses, metallicity and distance while the
star’s age has proven difficult to constrain as its colour and magnitude
change little during the main sequence phase.
The models outlined above don’t take extinction into account and the following
tests on simulated data assume no extinction. However, when applying the
method on real data, extinction can have a significant effect on a star’s
colour and magnitude such that, if left unchecked, parameter inference will be
highly inaccurate. To overcome this, extinction curves for each band are
applied. These are calculated according to a polynomial function of the star’s
colour and extinction at 550 nm which is provided for each star in the Gaia
catalogue using models from Fitzpatrick (1999). After calculating the
extinction in each band, this is simply subtracted from the measured
magnitude. This new magnitude is then used in SBI to recover astrophysical
parameters and extinction can be safely ignored.

##  3 Testing of the Model

###  3.1 Initial Test with Simulated Data

In order to test how this method is performing overall, SBI was run on a
sample of 1 million simulated stars with all parameters sampled from uniform
distributions. Each source was sampled 2,000 times and the median values of
these samples were compared to the simulated ‘true’ values. One way of
comparing the simulated parameter distributions is with a probability-
probability (P-P) plot. This plots the normalised cumulative distribution of
each recovered parameter against the true cumulative distribution (Gibbons &
Chakraborti, 2014). If the distributions agree perfectly, this plot follows a
straight 1-1 line. A P-P plot for the primary mass, mass ratio, age,
metallicity and distance distributions is shown inFigure 5.

Figure 5: P-P plot of the cumulative distributions of the recovered values
against true values. Dashed line is 1-1 line where the distributions are the
same.

The P-P plot demonstrates that the recovered distributions of primary mass and
distance correspond almost exactly with the simulated distribution. By
contrast, the recovered age distribution is very different as this parameter
is the most difficult to constrain by photometric methods. The P-P plot for
the mass ratio is almost always above the 1-1 plot which indicates this
parameter is frequently underestimated. The metallicity plot is mostly below
the 1-1 plot indicating this parameter is usually overestimated. This could be
due to the degeneracies when mapping parameters to observed magnitudes. There
could, therefore, be low metallicity binaries which are incorrectly labelled
as high metallicity single stars. The distributions of each parameter are
shown in detail in Figure 6 which plots 2D histrograms of recovered medians
against simulated values.

Figure 6: Comparison histograms between the simulated values of the parameters
and recovered medians. Density is on a log scale. The dashed line represents
1-1 mapping between simulated and recovered values.

As shown in Figure 6, similar to the P-P plot, strong correlation can be seen
between true and recovered values of primary mass and distance. Most of the
recovered mass ratio values are close to the true values but more than half of
them have been underestimated. Most of the recovered values of metallicity
agree with the simulated values but there are many for which the metallicity
has been overestimated. As stated above, this is likely due to degeneracies
and could also be due to the fact that the training set is denser at higher
metallicities leading the SBI recovery to favour stars at these higher
metallicities.

###  3.2 Recovered Binary Fractions

Of particular scientific interest is the prevalence of binary systems in
stellar populations. For the purposes of this study, a ‘binary’ is defined as
a source with $q\geq 0.5$. Below this mass ratio, it is unlikely a star’s
binarity can be determined by photometry as the magnitudes are not
sufficiently brighter than the single star case (as indicated in Figure 1.) A
‘true’ binary is defined as a simulated source with true $q\geq 0.5$ and a
‘recovered’ binary as a source with recovered median $q\geq 0.5$. Using this
definition of a binary, the relationship between recovered binary fraction and
colour and magnitude was investigated. From this, it is possible to identify
the regions of the H-R diagram at which this method is most sensitive and
accurate as well as those regions in which single stars are incorrectly
recovered as binaries. The binary fractions for the sample of true singles
(true $q<0.5$) and true binaries are shown on H-R diagrams in Figure 7.

(a) Binary fraction for stars with simulated $q<0.5$

(b) Binary fraction for stars with simulated $q>0.5$

Figure 7: Fractions of simulated single stars and binary stars recovered as
binaries shown on the H-R diagram.

Figure 7(b) shows a clear ‘stripe’ along the main sequence along which
binaries are successfully recovered. This demonstrates the limitation of
photometry as, below this stripe, the secondary star is not bright enough to
significantly affect the colour and magnitude. The stripe also gets less
pronounced at the lower mass end of the main sequence, indicating fewer
binaries are recovered for primary stars at these masses. This is probably due
to companions not significantly affecting the magnitude until $q$ is very high
(>0.8) and the larger magnitude uncertainties associated with fainter stars.
It should be noted that the binary fraction diagram for true singles shown in
Figure 7(a) also exhibits a stripe along the main sequence, however it is
confined to brighter magnitudes and is less pronounced. This stripe is due to
degeneracies in the relation between parameters and magnitudes. Many of the
single stars falsely recovered as binaries have high mass and/or metallicity
which, when an uncertainty is added, can make the star ‘look’ like a binary.
Both H-R diagrams show a region at the top end of the main sequence (and where
stars are moving towards the giant branch) in which almost 100% of stars are
recovered as binaries. In this region, there is very little difference between
Figures 7(a) and 7(b). This is an example of where the method fails, which is
most likely caused by the relative sparsity of sources in the training set as
shown in Figure 2.

###  3.3 Non-Single Star Catalogue

In order to further test the sensitivity of the recovery of binary fraction,
the SBI method was also run on the non-single star tables (NSS) from Gaia DR3.
These are stars which have been identified as astrometric, spectroscopic and
eclipsing binaries. For each star in the catalogue, the SBI was run as
previously described and calculated which sources had median recovered $q\geq
0.5$. The binary fraction on an H-R diagram is shown in Figure 8.

Figure 8: Binary fraction on an H-R diagram for Gaia’s NSS catalogue.

The result from the non-single star catalogue reveals a region of the H-R
diagram where binaries differ from singles so little in their photometry that
they are undetectable. Most of these would have low mass ratios or are at the
fainter end of the main sequence where magnitude uncertainties are larger. All
results from Figures 7 and 8 show a region at the brighter end of the main
sequence where all stars are recovered as binaries. From these results, a
region of the H-R diagram can then be identified in which the method is most
effective. When running the method on a large sample of Gaia data, this study
can then be restricted to this region in which the method sensitive to
binaries with minimal contamination from ‘false’ binaries.

##  4 Large Gaia Sample

When running SBI on a large sample of Gaia data, the sample is restricted to a
stripe along the main sequence in which this method has shown to be effective
at recovering binary fractions. This region is defined by a parallelogram on
the H-R diagram restricted to stars with Bp-Rp between 0–1.5 with absolute G
magnitude between two lines 2.2 magnitudes apart with a slope of 4.5. The
region studied is shown in Figure 9.

Figure 9: Region of the H-R diagram in which this study is focused, shown as a
red-dashed box on the distribution of the training set shown in Figure 2.

Additionally, conditions are set on the magnitude and distance, only
considering sources with apparent G$<$20 and distance$<$1,000 pc.

###  4.1 SBI Analysis and Binary fractions

The SBI method was run on Gaia sources in this region across the entire sky.
The sample of Gaia sources was distributed evenly across the sky, making it
possible to investigate whether the binary fraction changes as a function of
spatial position. This is shown in galactic coordinates in Figure 10.

Figure 10: Colour map of binary Fraction across the entire sky. This is
presented in galactic coordinates centred at $(l,b)=(0,0)$ with $l$ increasing
to the left. There is a clear region near the galactic centre where the binary
fraction is maximised and a bright point near $l=170^{\circ},b=-15^{\circ}$
which corresponds to the Taurus Molecular Cloud.

The binary fraction map in Figure 10 shows a concentration of binaries toward
the galactic centre and a lack of binaries in the galactic plane near
$l=180^{\circ}$. The perceived overabundance of binaries toward the galactic
centre could be due to the $\rho$-Ophiuchus star forming region and correlates
closely with the CO map from Dame & Thaddeus (2022). For sources at larger
distances ($>$300 pc), the overabundance can be explained by Gaia selection
effects. Gaia is likely to select only the brightest stars in this dense
region, which are more likely to be binaries. Outside the galactic plane,
there is little variation in binary fraction. However, there is a local peak
in binary fraction near galactic coordinates $l=170^{\circ},b=-15^{\circ}$ on
the far left of Figure 10. These coordinates correspond to the location of the
Taurus Molecular Cloud which, along with $\rho$-Ophiuchus, is one of the
closest star forming regions. These results suggest higher binary fraction in
close stellar associations.
The overall binary fraction is also expected to vary as a function of
metallicity (Moe et al., 2019). Figure 11 shows a 2D histogram comparing the
recovered values of mass ratio and metallicity as well as the binary fraction
as a function of metallicity.

(a) Histogram of recovered $q$ and [Fe/H]

(b) Binary Fraction as a Function of Metallicity

Figure 11: Histogram of mass ratio and metallicity and binary fraction as a
function of metallicity. For higher metallicities, result is strongly anti-
correlated.

The histogram in Figure 11(a) demonstrates the fact that most of the recovered
metallicities are near solar matallicities and values below $\sim$-0.5 are
rare. The mass ratio for high metallicity values peaks at low $q$ ($\sim$0.1)
and for matallicities $\sim$-0.5 it peaks at $q>0.6$. This is reflected in the
binary fraction plot shown in Figure 11(b) which shows anti-correlation
between binary fraction and metallicity which has been identified by previous
studies (Moe et al., 2019). However, below [Fe/H]$\sim$-0.6, the binary
fraction appears to increase with increasing metallicity and it is determined
that binaries with [Fe/H]<-1 are quite rare. However, as noted in Section 3.1,
some of the recovered single stars at medium metallicity ([Fe/H]$\sim$-1)
could actually be binary stars at lower metallicity. These incorrectly
characterized stars account for less then 5% of the total sample but could
significantly bias the result for low metallicities in which are relatively
uncommon. Overall, the binary fraction in the local stellar neighbourhood was
calculated to be 23.26$\pm$0.02%.

###  4.2 Combination with Paired Catalogue

As an additional study with this data set, the sources were cross-matched with
those in the paired catalogue with accurate radial velocity measurements
(Chance et al., 2022). The results of the SBI analysis detailed in Section 4.1
were combined with radial velocity amplitudes from paired to produce a period
distribution. Assuming these systems are on circular orbits, we can make use
of the relation:

$P=\frac{2\pi G}{(M_{1}+M_{2})^{2}}\left(\frac{M_{2}\rm{sin}i}{K}\right)^{3},$
(5)

where $M_{1}$ and $M_{2}$ are the primary and secondary masses of the system,
represented as distributions from the SBI analysis where $M_{2}=qM_{1}$. The
radial velocity $K$ is sampled from a Gaussian distribution of the same size
with a mean equal to the median $K$ value from Chance et al. (2022) and
standard deviation taken from the 16th and 84th percentiles of $K$. The
inclination angle $i$ is unknown, so a uniform prior distribution in $\cos{i}$
is adopted. The distributions of the recovered mass ratio and calculated
periods are shown in Figure 12(a). At low mass ratio ($q<0.2$) the calculation
of the period distribution is skewed towards extremely low values (as expected
by Equation 5) are likely unphysical. For $q>0.2$, the period distribution
maintains a constant shape which is shown in Figure 12(b).

(a) Distributions of recovered $q$ and calculated $P$

(b) Distribution of calculated $P$ for $q>0.2$.

Figure 12: Comparison of distributions of mass ratio and period distribution
excluding low $q$ values.

The period distribution shown in Figure 12(b) peaks at periods of slightly
more than 102 days which differs from previous studies outlined in Moe & Di
Stefano (2017) which have the distribution peaking at over $10^{4}$ days. This
is most likely due to the sensitivity of the radial velocity survey shown in
Figure 2 of Chance et al. (2022) which demonstrates this radial velocity
survey is not sensitive to companions beyond $\sim$10 AU.

##  5 Study of the Star Cluster M67

The binary fraction of star clusters is of particular interest as this can
provide information about stellar interactions and evolution. Additionally,
clusters that have been well-studied will already have constraints on the
astrophysical properties of its members. These constraints can then be applied
to the SBI results to improve their accuracy. In this section, an analysis of
367 stars in the open cluster M67 is presented. At an age of $\sim$4 Gyr, this
is one of the oldest known open clusters. With many Sun-like stars, it is also
one of the most well-studied clusters.

###  5.1 Forced Metallicity Distribution

This SBI method was run on a selection of Gaia sources in M67, identified by
position, proper motion and parallax measurements in the ranges specified by
Ghosh et al. (2022). The recovery was run as detailed in Section 2 which
produced a set of distributions for each star similar to that shown in Figure
4. After the properties of each star have been recovered in this way, a
metallicity distribution inspired by Önehag et al. (2014) is applied in which
a mean [Fe/H] of 0.06 and standard deviation of 0.1 are assumed. For each
star, every sample in the SBI recovery with a metallicity outside this
distribution was removed. An comparison for an example in the cluster is shown
in Figure 13 which plots the original result from SBI and the resultant
distributions when applying this forced metallicity distribution.

Figure 13: Comparison between parameter distributions before and after
applying an assumed metallicity distribution for the cluster. Shown for one
example in M67.

The new distributions shown in Figure 13 demonstrate a source which previously
had a broad mass ratio distribution that has now been more tightly constrained
an centered on slightly lower $q$ than before. This is due to the fact that,
for this particular star, the SBI recovery slightly favoured lower
metallicities and higher mass ratios which would produce similar magnitudes.
By forcing the metallicity to conform to what was already known about the
cluster, some of the degeneracy has been eliminated.

###  5.2 Mass ratios and Binary Fractions

This method of a forced metallicity distribution was applied to all stars in
the sample, producing more constrained estimates on stellar masses and
allowing a more accurate estimation of the binary fraction. An H-R diagram
(with extinction removed) showing the median mass ratios for stars in the
cluster is shown in Figure 14. The distinction between single and binary stars
is clearly evident both in the stars’ positions on the H-R diagram and the
recovered median masses.

Figure 14: H-R diagram showing median mass ratio from M67. There is a clear
distinction between the single star main sequence and the set of binaries
slightly above.

Using the previous definition of a binary (median $q$>0.5), the overall binary
fraction of this sample is found to be 34$\pm$3%. This is significantly higher
than the binary fraction calculated for the local stellar neighbourhood in
Section 4.1. The relation between this measured binary fraction and position
within the cluster was also investigated. This is represented in Figure 15
which shows the binary fraction as a function of angular separation from the
centre of the cluster. The centre is defined here as
(R.A.,Dec)=(132.8∘,+11.75∘.

Figure 15: Binary fraction as a function of angular distance from the centre
of the cluster. The binary fraction appears to follow a downward trend
indicating binaries are more common in the cluster’s interior but the large
uncertainties prevent any definite conclusion.

Plotting the binary fraction as a function of separation indicates a slight
downward trend which would be expected, based on previous studies indicating
binaries are more common in the interior of the cluster. However, due to the
relatively low number of sources in each bin, the results have high
uncertainties and the differences between binary fractions are not significant
enough for us to say this trend is real. Using the measured core radius of
8.24’ (Davenport & Sandquist, 2010) and the same definition of the cluster
centre, it was possible to measure the overall binary fraction inside and
outside the core. In total, 41 binaries were identified in the core out of a
total sample size of 92 giving a binary fraction of 45$\pm$8%. Outside the
core, 82 binaries were identified from a sample size of 275 giving a binary
fraction of 30$\pm$4%. While a higher binary fraction was found towards the
centre of the cluster, the large uncertainties on both binary fractions again
prevent a definite conclusion. It is possible the higher binary fraction near
the centre is due to a Gaia selection effect.

##  6 Summary and Conclusions

This work has shown that it is possible to recover parameters of main sequence
stars through photometry alone. The primary mass and metallicity can be
recovered to a high degree of accuracy while the age has proven difficult to
constrain. The binary mass ratio can also be recovered accurately in most
cases and, when applied to real data, can be used to search for binaries.
Using an assumed definition of a binary, it is then possible to estimate the
binary fraction of a sample. Using a simulated data set, regions in
colour/magnitude space were identified in which this method produces more
accurate results. This is possibly due to the higher density of the training
set in these regions.
Focusing in on the region on the H-R diagram where the method produces
accurate results, it was possible to analyse the binary fractions of the local
stellar neighbourhood as functions of position on the sky and metallicity.
Overall, the binary fraction was found to change little with position in the
sky but this study has reproduced the expected anti-correlation between binary
fraction and metallicity. It has been shown that these results can be combined
with radial velocity studies to produce period distributions with a shape that
agrees with previous studies but peaks at lower periods than expected.
However, this can be attributed to lack of sensitivity at wide separations.
Through the study of M67, it has been shown that applying a known distribution
of a star’s metallicity can greatly improve the precision of the recovery of
other parameters including mass ratio. A binary fraction was calculated which
is almost constant regardless of distance from the centre but exhibits a
slight downward trend. Out binary fraction for M67 was calculated to be
34$\pm$3% which is higher than the 23.26$\pm$0.02% calculated for field stars.
This suggested a similar result to previous studies (Offner et al., 2022)
which found that binaries are more common in environments with higher stellar
density.
Moving forward, this project enables a more complete picture of the binary
fraction of sources in the Gaia data. Radial velocity and astrometric surveys
are primarily concerned with detecting dark companions and, with this extra
knowledge of unreolved bright companions, these surveys will be able to place
better constraints on hidden companions.

## Data Availability

The code used to recover astrophysical parameters with SBI is available in a
github repository (https://github.com/awallace142857/sbi_code) which also
includes some Gaia data. Due to the size of the data files, the remaining data
are available from the corresponding author on reasonable request.

## Acknowledgements

A.L.W. acknowledges the contributions of Andrew Casey who helped start this
project and edit this paper in its early stages. A.L.W. also thanks Anthony
Brown, Daniel Foreman-Mackey and Quadry Chance for their helpful suggestions
regarding the discussion and overall science goals of the study. This work was
funded by the Australian Research Council Discovery Grant DP210100018.

## References

  * Badenes et al. (2018) Badenes C., et al., 2018, ApJ, 854, 147
  * Chance et al. (2022) Chance Q., Foreman-Mackey D., Ballard S., Casey A. R., David T. J., Price-Whelan A. M., 2022, paired: A Statistical Framework for Determining Stellar Binarity with Gaia RVs. I. Planet Hosting Binaries, doi:10.5281/zenodo.6765903, https://doi.org/10.5281/zenodo.6765903
  * Clark et al. (2012) Clark B. M., Blake C. H., Knapp G. R., 2012, ApJ, 744, 119
  * Cranmer et al. (2020) Cranmer K., Brehmer J., Louppe G., 2020, Proceedings of the National Academy of Sciences, 117, 30055
  * Dame & Thaddeus (2022) Dame T. M., Thaddeus P., 2022, ApJS, 262, 5
  * Davenport & Sandquist (2010) Davenport J. R. A., Sandquist E. L., 2010, The Astrophysical Journal, 711, 559
  * Deacon & Kraus (2020) Deacon N. R., Kraus A. L., 2020, MNRAS, 496, 5176
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, Annual Review of Astronomy and Astrophysics, 51, 269
  * Duchêne et al. (2018) Duchêne G., Lacour S., Moraux E., Goodwin S., Bouvier J., 2018, Monthly Notices of the Royal Astronomical Society, 478, 1825
  * Fitzpatrick (1999) Fitzpatrick E. L., 1999, PASP, 111, 63
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Ghosh et al. (2022) Ghosh E. M., Sulistiyowati Tucio P., Fajrin M., 2022, Journal of Physics: Conference Series, 2214, 012009
  * Gibbons & Chakraborti (2014) Gibbons J., Chakraborti S., 2014, Nonparametric Statistical Inference, Fourth Edition: Revised and Expanded. Taylor & Francis, https://books.google.com.au/books?id=kJbVO2G6VicC
  * Hahn & Melchior (2022) Hahn C., Melchior P., 2022, The Astrophysical Journal, 938, 11
  * Hartigan & Kenyon (2001) Hartigan P., Kenyon S., 2001, Symposium - International Astronomical Union, 200, 496–500
  * Lindegren et al. (2021) Lindegren L., et al., 2021, Astronomy & Astrophysics, 649, A4
  * Mazzola et al. (2020) Mazzola C. N., et al., 2020, Monthly Notices of the Royal Astronomical Society, 499, 1607
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Moe et al. (2019) Moe M., Kratter K. M., Badenes C., 2019, The Astrophysical Journal, 875, 61
  * Morton (2015) Morton T. D., 2015, isochrones: Stellar model grid package, Astrophysics Source Code Library, record ascl:1503.010 (ascl:1503.010)
  * Mugrauer et al. (2022) Mugrauer M., Zander J., Michel K.-U., 2022, Astronomische Nachrichten, 343, e24017
  * Offner et al. (2022) Offner S. S. R., Moe M., Kratter K. M., Sadavoy S. I., Jensen E. L. N., Tobin J. J., 2022, arXiv e-prints, p. arXiv:2203.10066
  * Saleh & Rasio (2009) Saleh L. A., Rasio F. A., 2009, ApJ, 694, 1566
  * Traven et al. (2020) Traven G., et al., 2020, Astronomy & Astrophysics, 638, A145
  * Önehag et al. (2014) Önehag A., Gustafsson B., Korn A., 2014, Astronomy & Astrophysics, 562


#  ChatGPT, Large Language Technologies, and the Bumpy Road of Benefiting
Humanity

Atoosa Kasirzadeh In DailyNous — March 16, 2023 — “Philosophers on Next-
Generation Large Language Models” series, edited by Annette Zimmermann.
Contributions by Abeba Birhane, Atoosa Kasirzadeh, Finton Mallory, Regina
Rini, Eric Schwitzgebel, Luke Stark, Karina Vold, and Annette Zimmermann.

From tech moguls in Silicon Valley to those who have the luxury of indulging
in the exploration of cutting-edge AI technologies, OpenAI’s ChatGPT has
captured the imagination of many with its conversational AI capabilities. The
large language models that underpin ChatGPT and similar language technologies
rely on vast amounts of textual data and alignment procedures to generate
responses that can sometimes leave users pondering whether they’re interacting
with a piece of technology or a human. While some view making language agents
such as ChatGPT merely as a significant step in developing AI for linguistic
tasks, others view it as a vital milestone in the ambitious pursuit of
achieving artificial general intelligence – AI systems that are generally more
intelligent than humans. In a recent blogpost OpenAI’s CEO, Sam Altman,
emphasizes the ambitious role of this technology as a step towards building
“artificial general intelligence” that “benefits all of humanity.”

ChatGPT promises to enhance efficiency and productivity with its remarkable
capabilities. One impressive feature is its ability to summarize texts. For
example, if you do not have time to read Sam Altman’s complex argument from
2018 when he agreed with Emily Bender, a prominent linguist from the
University of Washington, that humans are not stochastic parrots, you can ask
ChatGPT and it will summarize the argument in a blink of an eye:

Or if you are curious to have a summary of David Chalmers’ 2019 speech at the
United Nations about the dangers of virtual reality, ChatGPT comes to your
service:

Impressive outputs, ChatGPT! For some people, these results might look like
watching a magician pull a rabbit out of a hat. However, we must address a few
small problems with these two summaries: the events described did not happen.
Sam Altman did not agree with Emily Bender in 2018 about humans being
stochastic parrots; the discussion regarding the relationship between
stochastic parrots, language models, and human’s natural language processing
capacities only got off the ground in a 2021 paper “on the dangers of
stochastic parrots: can language models be too big?” [1]. Indeed, in 2022
Altman tweeted that we are stochastic parrots (perhaps sarcastically).

Similarly, there is no public record of David Chalmers giving a speech at the
United Nations in 2019. Additionally, the first arXiv link in the bibliography
takes us to the following preprint, which is neither written by David Chalmers
nor is titled “The Dangers of Stochastic Parrots: Can Language Models Be Too
Big?”

The second bibliography link takes us to a page that cannot be found.

These examples illustrate that outputs from ChatGPT and other similar language
models can include content that deviates from reality and can be considered
hallucinatory. While some researchers may find value in the generation of such
content, citing the fact that humans also produce imaginative content, others
may associate this with the ability of large language models to engage in
counterfactual reasoning. However, it is important to recognize that the
inaccuracies and tendency of ChatGPT to produce hallucinatory content can have
severe negative consequences, both epistemically and socially. Therefore, we
should remain cautious in justifying the value of such content and consider
the potential harms that may arise from its use.

One major harm is the widespread dissemination of misinformation and
disinformation, which can be used to propagate deceptive content and
conspiracies on social media and other digital platforms. Such misleading
information can lead people to hold incorrect beliefs, develop a distorted
worldview, and make judgments or decisions based on false premises. Moreover,
excessive reliance on ChatGPT-style technologies may hinder critical thinking
skills, reduce useful cognitive abilities, and erode personal autonomy. Such
language technologies can even undermine productivity by necessitating
additional time to verify information obtained from conversational systems.

I shared these two examples to emphasize the importance of guarding against
the optimism bias and excessive optimism regarding the development of ChatGPT
and related language technologies. While these technologies have shown
impressive progress in natural language processing, their uncontrolled
proliferation may pose a threat to the social and political values we hold
dear.

I must acknowledge that I am aware and excited about _some_ potential benefits
of ChatGPT and similar technologies. I have used it to write simple Python
codes, get inspiration for buying unusual gifts for my parents, and crafting
emails. In short, ChatGPT can undoubtedly enhance some dimensions of our
productivity. Ongoing research in AI ethics and safety is progressing to
minimize the potential harms of ChatGPT-style technologies and implement
mitigation strategies to ensure safe systems.111For two examples, see Taxonomy
of Risks posed by Language Models [3] for our recent review of such efforts as
well as Anthropic’s Core Views on AI safety. These are all promising
developments.

However, despite _some_ progress being made in AI safety and ethics, we should
avoid oversimplifying the promises of artificial intelligence “benefiting all
of humanity”. The alignment of ChatGPT and other (advanced) AI systems with
human values faces numerous challenges.222For a philosophical discussion, see
our paper In conversation with Artificial Intelligence: aligning language
models with human values [2]. One challenge is that human values can be in
conflict with one another. For example, we might not be able to make
conversational agents that are simultaneously maximally helpful and maximally
harmless. Choices are made about how to trade-off between these conflicting
values, and there are many ways to aggregate the diverse perspectives of
choice makers. Therefore, it is important to carefully consider which values
and whose values we align language technologies with and on what legitimate
grounds these values are preferred over other alternatives.

Another challenge is that while recent advances in AI research may bring us
closer to achieving some dimensions of human-level intelligence, we must
remember that intelligence is a multidimensional concept. While we have made
great strides in natural language processing and image recognition, we are
still far from developing technologies that embody unique qualities that make
us human–our capacity to resist, to gradually change, to be courageous, and to
achieve things through years of dedicated effort and lived experience.

The allure of emerging AI technologies is undoubtedly thrilling. However, the
promise that AI technologies will benefit all of humanity is empty so long as
we lack a nuanced understanding of what humanity is supposed to be in the face
of widening global inequality and pressing existential threats. Going forward,
it is crucial to invest into rigorous and collaborative AI safety and ethics
research. We also need to develop standards in a sustainable and equitable way
that differentiate between merely speculative and well-researched questions.
Only the latter enable us to co-construct and deploy the values that are
necessary for creating beneficial AI. Failure to do so could result in a
future, where our AI technological advancements outstrip our ability to
navigate the ethical and social implications. This path we do not want to go
down.

## References

  * Bender et al. [2021] Bender, E. M., Gebru, T., McMillan-Major, A., and Shmitchell, S. (2021).  On the Dangers of Stochastic Parrots: Can Language Models Be Too Big?  In _Proceedings of the 2021 ACM Conference on Fairness, Accountability, and Transparency_. 610–623.
  * Kasirzadeh and Gabriel [2022] Kasirzadeh, A. and Gabriel, I. (2022).  In conversation with Artificial Intelligence: aligning language models with human values.  _arXiv preprint arXiv:2209.00731_ .
  * Weidinger et al. [2022] Weidinger, L., Uesato, J., Rauh, M., Griffin, C., Huang, P.-S., Mellor, J., Glaese, A., Cheng, M., Balle, B., Kasirzadeh, A. et al. (2022).  Taxonomy of risks posed by language models.  In _2022 ACM Conference on Fairness, Accountability, and Transparency_. 214–229.


# Knowledge-Assisted Reasoning
of Model-Augmented System Requirements with Event Calculus and Goal-Directed
Answer Set Programming111UT Dallas authors are supported by NSF (grants IIS
1718945, IIS 1910131, IIP 1916206), Amazon and DoD. Support for this project
also comes from the FIT BUT internal project FIT-S-20-6427 and by the H2020
ECSEL project Arrowhead Tools. We would also like to dedicate our work to the
memory of second author’s Father.

Brendan Hall1 1Honeywell Advanced Technology, Plymouth, USA2The University of
Texas at Dallas, Richardson, USA3Honeywell International s.r.o & Brno Univ. of
Technology, Brno, Czech Republic4Universidad Rey Juan Carlos, Madrid, Spain
Sarat Chandra Varanasi2 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Jan Fiedor3 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Joaquín Arias4 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Kinjal Basu2 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Fang Li2 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Devesh Bhatt1 1Honeywell Advanced Technology, Plymouth, USA2The
University of Texas at Dallas, Richardson, USA3Honeywell International s.r.o &
Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan Carlos,
Madrid, Spain Kevin Driscoll1 1Honeywell Advanced Technology, Plymouth,
USA2The University of Texas at Dallas, Richardson, USA3Honeywell International
s.r.o & Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan
Carlos, Madrid, Spain Elmer Salazar2 1Honeywell Advanced Technology, Plymouth,
USA2The University of Texas at Dallas, Richardson, USA3Honeywell International
s.r.o & Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan
Carlos, Madrid, Spain Gopal Gupta2 1Honeywell Advanced Technology, Plymouth,
USA2The University of Texas at Dallas, Richardson, USA3Honeywell International
s.r.o & Brno Univ. of Technology, Brno, Czech Republic4Universidad Rey Juan
Carlos, Madrid, Spain 1,3{brendan.hall, jan.fiedor, devesh.bhatt,
<EMAIL_ADDRESS>2{sarat-chandra.varanasi, kinjal.basu, fang.li,
ees101020<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider requirements for cyber-physical systems represented in constrained
natural language. We present novel automated techniques for aiding in the
development of these requirements so that they are consistent and can
withstand perceived failures. We show how cyber-physical systems’ requirements
can be modeled using the event calculus (EC), a formalism used in AI for
representing actions and change. We also show how answer set programming (ASP)
and its query-driven implementation s(CASP) can be used to directly realize
the event calculus model of the requirements. This event calculus model can be
used to automatically validate the requirements. Since ASP is an expressive
knowledge representation language, it can also be used to represent contextual
knowledge about cyber-physical systems, which, in turn, can be used to find
gaps in their requirements specifications. We illustrate our approach through
an altitude alerting system from the avionics domain.

##  1 Introduction

Developing effective requirements is crucial for success in building a system.
The earlier the requirements are validated, the fewer problems system
developers will encounter later. Current automation of requirements
engineering tasks attempt to ensure their consistency and adequacy, namely,
that they can withstand perceived failures. However, such automated support
remains limited. In this work, we present novel automated techniques for
aiding the development of model-augmented requirements that are adequate—to
the extent that adequacy can be established—and consistent. Thus, we can have
more confidence in the requirements. We limit ourselves to requirements for
cyber-physical systems, particularly those in avionics. We assume that
requirements are generated within the MIDAS (Model-Assisted Decomposition and
Specification) [12] environment, and are expressed in CLEAR (Constrained
Language Enhanced Approach to Requirements) [11], a constraint natural
requirement language based on EARS [15].

Our main contribution in this work is to show how the Event Calculus (EC) [20,
21] and Answer Set Programming (ASP) [10] can be used to formalize constrained
natural language requirements for cyber-physical systems and perform
knowledge-assisted reasoning over them. ASP is a logic-based knowledge
representation language that has been prominently used in AI. Our work builds
upon recent advances made within the s(CASP) system [3], a query-driven (or
goal-directed) implementation of predicate ASP that supports constraint
solving over reals, permitting the faithful representation of time as a
continuous quantity. The s(CASP) system permits the modeling of event calculus
directly [4]. An advantage of using the event calculus—in contrast to automata
and Kripke structure-based approaches—is that it can directly model cyber-
physical systems, thereby avoiding “pollution" due to (often premature) design
decisions that must be made otherwise. The event calculus is a formalism—a set
of axioms—for modeling dynamic systems and was proposed by artificial
intelligence researchers to solve the frame problem [20, 21]. Our methods have
been developed within the MIDAS framework of Honeywell [12]. The MIDAS
framework embodies the essence of the Property Management Model [17] within a
model-augmented requirement ecosystem, which integrates constrained natural
language requirements within a layered information model expressed in Object
Process Model (OPM) [8]. The details of MIDAS, PMM and OPM are not important
here, instead, the primary goal of this work is to explore how constrained
natural language requirements, specified within MIDAS, using the CLEAR
notation, can be automatically reasoned about and analyzed using the event
calculus and query-driven answer set programming. Specifically, we explore:

  1. 1.

How ASP-based model checking (over dense time) can validate specified system
behaviors wrt system properties.

  2. 2.

How application of abductive reasoning can extend ASP-based model checking to
incorporate domain knowledge and real-world/environmental
assumptions/concerns.

  3. 3.

How knowledge-driven analysis can identify typical requirement specification
errors, and/or requirement constructs which exhibit areas of
potential/probable risk.

The rest of the paper is organized as follows. Section 2 introduces the
enabling technologies of answer set programming (ASP), a goal-directed
implementation of ASP called s(CASP), and the event calculus. Section 3
discusses the CLEAR notation for specifying requirements, and illustrates it
with an example. Section 4 presents an application developed within Honeywell
Corporation using our approach. We also illustrate requirement defect
discovery using s(CASP) for property-based model-checking as well as discuss
how more general knowledge of potential requirements defects may detect
defects that traditional techniques may not be able to find. Section 5
presents our conclusions and future work.

The main contribution of our research is to show how the event calculus and
goal-directed answer set programming together serve as a promising framework
for modeling and reasoning over cyber-physical and avionics systems’
requirements. We also demonstrate the utility of abductive reasoning to find
the assumptions under which a property holds or does not hold (along with
automatically providing justification in both cases) for requirements
specification. Given that ASP is a knowledge representation language, we also
illustrate how knowledge about system context (e.g., processing platform may
malfunction causing them to be reset) can be utilized to automatically find
gaps in requirements as well.

##  2 Enabling Technologies

###  2.1 Answer Set Programming

Answer Set Programming (ASP) is a declarative paradigm that provides the
stable model semantics for logic programs (with negation-as-failure). ASP is a
highly expressive paradigm that can elegantly express complex reasoning
methods, including those used by humans, such as default reasoning, deductive
and abductive reasoning, counterfactual reasoning, constraint satisfaction,
etc. [5, 10].

The rules in an ASP program are of the form: p :- q1, ..., qm, not r1, ...,
not rn. where $m\geq 0$ and $n\geq 0$. Each of p and qi ($\forall i\leq m$) is
a literal, and each not rj ($\forall j\leq n$) is a naf-literal (not is a
logical connective called negation-as-failure or default negation). The
literal not rj is true if proof of rj fails. Negation as failure allows us to
take actions that are predicated on failure of a proof. Thus, the rule r :-
not s. states that r can be inferred if we fail to prove s. Note that in the
rule above, p is optional. Such a headless rule is called a constraint, which
states that conjunction of qi’s and not rj’s should yield false. Thus, the
constraint :- u, v. states that u and v cannot be both true simultaneously in
any model of the program (a model is called an answer set).

ASP also allows for assumption-based reasoning via abduction [13]. Abduction
is a form of reasoning concerned with the generation and evaluation of
explanatory hypotheses. In abductive reasoning, given the premise
$P\Rightarrow Q$, and the observation $Q$, one surmises or assumes (abduces)
that $P$ holds. More generally, given a theory $T$, an observation $O$, and a
set of abducibles $A$, then $E$ is an explanation of $O$ (where $E\subset A$)
if:

1\. $T\cup E\models O$

2\. $T\cup E$ is consistent

We can think of abducibles $A$ as a set of assumptions. Abduction allows us to
find “gaps" in our knowledge that must entail a property in order for that
property to hold.

###  2.2 The s(CASP) System

Considerable research has been done on answer set programming since the
inception of the stable model semantics that underlies it [10]. A major
problem with ASP implementations is that programs have to be grounded and SAT-
solver-based implementations such as CLASP [9] used to execute the
propositionalized program to find the answer sets. There are multiple problems
with this SAT-based implementation approach, which include exponential blowup
in program size, having to compute the entire model, and not being able to
produce a justification for a conclusion. Goal-directed implementations of ASP
[3], called s(ASP) and s(CASP), work directly on predicate ASP programs (i.e.,
no grounding is needed) and are query-driven (similar to Prolog). The s(ASP)
and s(CASP) systems only explore the parts of the knowledge-base that are
needed to answer the query, and they provide a proof tree that serves as
justification for the query. The s(ASP) and s(CASP) systems support predicates
with arbitrary terms as arguments as well as constructive negation [3].

In the work reported here, we will mainly use the s(CASP) system that
additionally supports constraint solving over reals, which is important for
reasoning faithfully about (continuous) time. The s(CASP) system is the key
technology for representing and analyzing CLEAR requirements modeled with the
event calculus. The s(CASP) system also directly supports abductive reasoning
and can provide justification for a query.

###  2.3 Event Calculus

The Event Calculus [20, 21] is a formalism for modeling dynamic systems which
originated from the field of AI where it was devised to solve the frame
problem. The event calculus is organized around events and fluents. We use the
event calculus formalism to model an evolving world in a cyber-physical
system. The EC models the world that is changing due to events that are
happening that, in turn, influence fluents, where a fluent is a variable that
changes with time such as location of an object or temperature of the boiler.
EC provides axioms for stating the conditions under which events happen or
conditions under which events initiate, terminate, or release a fluent. EC
also provides axioms for when a fluent holds. EC models the laws of inertia
that hold in the real world. For instance, the location (a fluent) of an
object does not change unless the event “move” happens. Event calculus uses
reasoning over time, as events are occurrences in time. The Event Calculus is
elegantly modeled using ASP, in particular, with the query-driven s(CASP)
system due to its support for real-time constraints.

The basic event calculus axioms are shown in Fig. 1. The happens/2 predicate
models the occurrence of an event $e$ at a particular time instance $t_{i}$,
while holdsAt/2 specifies the condition under which a fluent $f$ holds at time
$t_{j}$. The axioms are easy to follow. More details can be found elsewhere
[18, 4].

In a cyber-physical system, events correspond to actuator actions and fluents
correspond to effects of actuators. Fluents can also describe the conditions
under which events can occur. A fluent is an entity whose value may change
over time and, thus, in case of a CPS, they reflect changes in sensor readings
as well as internal component(s) state. Fluents may be impacted by events,
e.g, furnace shutoff (event) may lead to temperature (fluent) to drop. The
state of a system is represented by the values of all its fluents. The event
calculus simply models how events lead to changes in fluents. Fluents, in
turn, may trigger events, e.g., temperature reaching 80 degrees may cause
furnace shutoff. Thus, one can model a cyber physical system directly in the
event calculus without having to make any representational assumptions. This
is in contrast to other formalisms such as automata-based or Kripke structure-
based approaches where one has to first discern the states of the system
followed by the transitions between these states and then model it. The
discerning of states and transitions for these methods is done too early in
the design process leading to (possibly erroneous) assumptions being
incorporated into the design. In the event calculus, states are implicit, as a
state is defined as consisting of values of all fluents at any given moment.
Because of the implicit nature of system states, a state needs to be
constructed only at the time we wish to verify something, e.g., a property.

BEC1. $StoppedIn(t_{1},f,t_{2})\equiv$ $\exists\ e,t\ (Happens(e,t)\land\
t_{1}<t<t_{2}\land(Terminates(e,f,t)\lor Releases(e,f,t)))$
---
BEC2. $StartedIn(t_{1},f,t_{2})\equiv$ $\exists\ e,t(Happens(e,t)\land\
t_{1}<t<t_{2}\land(Initiates(e,f,t)\lor Releases(e,f,t)))$
BEC3. $HoldsAt(f_{2},t_{2})\leftarrow$ $Happens(e,t_{1})\land
Initiates(e,f_{1},t_{1})\land Trajectory(f_{1},t_{1},f_{2},t_{2})\land\lnot
StoppedIn(t_{1},f_{1},t_{2})$
BEC4. $HoldsAt(f,t)\leftarrow$ $InitiallyP(f)\land\lnot StoppedIn(0,f,t)$
BEC5. $\lnot HoldsAt(f,t)\leftarrow$ $InitiallyN(f)\land\lnot
StartedIn(0,f,t)$
BEC6. $HoldsAt(f,t_{2})\leftarrow$ $Happens(e,t_{1})\land
Initiates(e,f,t_{1})\land\ t_{1}<t_{2}\land\lnot StoppedIn(t_{1},f,t_{2})$
BEC7. $\lnot HoldsAt(f,t_{2})\leftarrow$ $Happens(e,t_{1})\land
Terminates(e,f,t_{1})\land\ t_{1}<t_{2}\land\lnot StartedIn(t_{1},f,t_{2})$
Figure 1: Formalization of BEC (Basic Event Calculus) axioms [18].

##  3 CLEAR Notation for Requirements

Constrained Language Enhanced Approach to Requirements (CLEAR) extends the
Easy Approach to Requirements Syntax (EARS) developed initially by Alistair
Mavin [15]. EARS is an industrially pragmatic approach based on using five
structured templates and keywords. Studies [16, 14] have shown that use of
EARS reduces requirements errors while improving requirement quality and
readability. The EARS keywords and templates are illustrated in Figure 2.

Although the EARS approach is often presented as a syntactic guide to writing
requirements, we have found the power of EARS to be more related to the
mindset than syntax. The different perspectives that are explored using EARS-
based keywords stimulate a lot of discussion relating to the core intent of
the requirement under consideration. The separation and clarification of event
and state semantics, using the WHILE and WHEN constructs are useful to
constrain and define the behavioral details of cyber-physical systems at
different levels of abstraction.

The EARS approach also calls for the separation of the nominal and off-nominal
specification and the IF keyword is reserved for defining the required
response to an unwanted or exceptional condition. In practice we have found
that it is often the insufficient understanding and naive assumptions
regarding system failure scenarios that lead to omissions of critical
requirements.

EARS designates optional and configurable behaviors/attributes with the WHERE
keyword. Once again, this is also an area of potential vulnerability regarding
requirement completeness. Significant hazards can emerge from
misconfiguration, as illustrated by the A400M incident in Seville [22]. Hence,
separating concerns helps address such issues in a more systematic fashion.

CLEAR embodies the full spirit of EARS and refines the terms used within EARS
constructs to improve specificity. CLEAR also introduces additional constructs
such as Upon Initialization to force a systematic examination of known problem
areas, such as system initialization. CLEAR is used today to support automated
analyses (e.g. requirement consistency checking) and software level test
generation [6]. Through MIDAS, it is intended to move and automate such
analysis to the system level.

###  3.1 Mapping Requirements to the Event Calculus

We discuss the mapping of CLEAR-based requirements [11] to the Extended Event
Calculus formalization of Shanahan[21]. CLEAR is a notation for expressing
requirements for cyber physical systems based on EARS (Easy Approach to
Requirement Syntax) widely adopted in the avionics industry. For simplicity,
we limit ourselves to discussing only EARS. EARS provides natural language
sentence templates for an engineer specifying a system to follow for writing
requirements. Keywords ‘WHEN’, ‘WHERE’ and ‘IF’-‘THEN’ are used in these
templates and play a major role. We make use of the extended event calculus
formalism because it enables events to have duration. For cyber-physical real-
time systems, response times are important, hence formally budgeting the
allocation of time throughout the levels of function & temporal decomposition
are primary concerns. Figure 2 summarizes the EARS notation. Details about
EARS can be found elsewhere [15, 14].

Figure 2: EARS Templates

Conceptually, the mapping of the extended event calculus to CLEAR is
relatively straightforward, given that EARS already informally separates event
and state semantics. To formalize the mapping, we only need to formally bound
the level of temporal abstraction over which the CLEAR requirements are
defined. The EARS keyword WHEN then needs to define the requirements in
relation to the upward or downward perspective, with the downward facing view
having the form of state invariant constraints, and the upward facing view
declaring actionable behaviours. For actionable behaviours, the level of
temporal abstraction also guides which of the EARS keywords to apply to the
specification as outlined below:

  * •

WHEN: discrete event-driven specification is used when a single observation at
the level of temporal abstraction is sufficient to verify the required
behavior. That is, WHEN establishes causal relationship.

  * •

WHILE: continuous state-driven specification is used when multiple
observations at the level of temporal abstraction are needed to verify the
required behavior. WHILE may also be used as a state qualifier for event-
driven specifications, in cases where there are stateful conditions guarding
the event.

Within MIDAS, functional influence manifests as the stateful changes to
objects that surround the functional intent boundary. This maps very cleanly
to the extended event calculus, once the state of the said objects at the
function boundary are encoded as fluents within the event calculus formalism.
Hence, to map the MIDAS CLEAR actionable requirements to the extended event
calculus, we simply perform the following (details of MIDAS are not important
here):

  1. 1.

Declare each MIDAS object state as a fluent.

  2. 2.

Declare happens predicate that characterize the discrete changes in state that
triggers the WHEN condition.

  3. 3.

Declare initiates and terminates predicates to bind the happens predicates to
the influenced external state changes, with holdsAt state qualifiers for all
WHILE preconditions .

  4. 4.

Declare trajectory predicates to characterize continuous changes that evolve
over multiple intervals of the temporal precision.

We illustrate the mapping of an EARS requirement related to aircraft design:

> While the aircraft is on-ground, when the requested door position becomes
> open, the door control system shall change the state of the cargo door from
> closed to open within 10 secs.

Mapping this to the EARS templates (Fig. 2) we see that the aircraft is on-
ground is a precondition, the requesting of the cargo door to open is the
trigger, and the opening of the cargo door within 10 seconds is the system
response, requiring the change of environmental state, with an assumed level
of temporal precision. In this requirement the functional intent manifests as
a change in the physical state of the cargo door, in response to a change in
state that reflects the pilot requesting the door to open. Hence, the event
calculus formalism of this requirement is as follows:

⬇

fluent(door_requested_position(RP)). fluent(door_state(DS)).

fluent(aircraft_state(AS)).

happens(pilot_requests_door_to_open, T) :-
-holdsAt(door_requested_position(open), T1),

holdsAt(door_requested_position(open), T2), T2 #=< T1+10, T1 #< T, T #=< T2.

initiates(pilot_requests_door_to_open, door_state(open), T) :-

holdsAt(aircraft_state(is_on_ground), T).
terminates(pilot_requests_door_to_open, door_state(closed),T) :-

holdsAt(aircraft_state(on_ground), T). door_response_is_correct_condition :-
-holdsAt(door_requested_position(open), TB), -holdsAt(door_state(open), TB),

holdsAt(aircraft_state(on_ground), TH), holdsAt(door_state(open), TE),

TB #=< TH, TB #< TE, TE #=< TB+10.

According to the semantics of the Extended Event Calculus [21], a fluent is
considered to be initiated at any time point during the event duration. This
is semantically consistent with the additional within constraint that the
CLEAR requirement notation [11] provides. Note also that the extended event
calculus formalism does not distinguish between the sensing of the
environmental change and the resulting time of establishing environmental
influence. This is also consistent with the level of functional and temporal
decomposition at the current level of abstraction. Hence, when validating and
discharging properties and constraints above this level of specification, the
event calculus model will explore the impact of the function influence
starting at any point within the event duration, which in our example is 10
seconds.

Note also that the initiated action, i.e., the changing of the door state, is
guarded by the aircraft state of ‘on ground’. Should this hold when the
respective triggering event happens, the initiates predicate will fail, and
therefore there will be no external functional influence. The mapping of IF is
the same as the WHEN, given that they are both conditioned on events. The
primary difference is that IF characterizes triggering off-nominal state
transitions such as component failure scenarios.

As noted above actionable WHILE requirements are used to characterize
continuous behaviours and/or constraints that need to be expressed over
multiple intervals at the current level of temporal abstraction. They map to
trajectories within the extended event calculus formalism. We show an example
and its mapping to EC below:-

> While the door is opening, the rate of change in door position shall
> increase positively with a maximum rate of 5 degrees per second.

⬇

trajectory(door_opening,T1,door_position(B),T2) :-

holdsAt(door_position(A),T1),

holdsAt(door_position(B),T2), A #< B, B #=< A+5, T2 #= T1+1.

door_rate_ok :- trajectory(door_opening,T1,door_position(P),T2),

holdsAt(door_opening,TH), T1 #=< TH, TH #=< T2.

ASP with EC can also be used to test the adequacy of the requirements to the
extent possible. In the case of avionics, the knowledge is simple and relates
to Single Event Upsets (or SEUs) [19] which can induce random resets of
avionic software systems. We can model the destructive power of the SEU using
a reset event. The reset event overrides the constructed internal state of the
system, forcing the initialization state to be re-established. Our goal in
introducing this reset is to move the designer to consider how robust the
system initialization logic is to such transient resets. In ASP, an extraneous
event (e.g., reset) is represented as an abducible (a reset may or may not
happen). We would want to know if our system will still behave correctly in
presence of this extraneous event. The knowledge assumes that the reset event
can only override the cyber system (software state) and does not affect the
continuous state of the physical world. Whether our system will behave
correctly or not can be established with the abductive reasoning supported by
s(CASP). This form of abductive reasoning has significant potential in helping
establish the completeness of requirements for cyber-physical systems.

##  4 Demonstration Case Study

To explore the feasibility of our approach, we applied it to several systems
under development at Honeywell. The primary goal of these case-studies is to
explore how specialist domain knowledge can be incorporated into the
requirement review process to help assess requirement sufficiency. Rather than
encode hard-coded knowledge of specific cases, we are striving to illustrate
the application of more general reusable knowledge. Our findings have been
encouraging, and we have been able to supply useful feedback to live
engineering teams. We illustrate our ideas with one of the systems developed.

###  4.1 Altitude Alerting System

The altitude alerting system should issue alerts about aircraft altitude error
thresholds in a timely fashion. Figure 3 shows the alerts that should be
turned on/off within one second as the absolute altitude crosses alert
thresholds. We first present the various notions associated with the altitude
alert system’s requirements, and show their mapping from CLEAR to Event
Calculus and ASP

Figure 3: Altitude Alerts w.r.t Altitude Error Thresholds

During the operation of the aircraft, the pilot selects a target altitude. The
alert system monitors the aircraft’s barometric altitude trajectory and issues
alerts when this barometric altitude becomes potentially hazardous in relation
to the selected altitude. If the error crosses certain threshold values, then
the system should promptly warn the pilot so that she can take appropriate
action. Similarly, the alerting system can also withdraw the alerts if other
threshold values are met. The two primary alerts are capture alert and
departure alert explained in Figure 3. In the following, lt and gt are
abbreviated for ‘less than’ and ‘greater than’, respectively. The predicate
error_becomes_lt(1000,T) asserts if the altitude error becomes less than 1000
ft. at time T.

REQ1:222It should be noted that these requirements are significantly simpler
than the original internal baseline which consisted of a state-based
specification. This simplification was enabled by MIDAS’ conceptual modeling
approach, which forces the specification to be expressed using conceptual
abstractions, over object states that surround the boundary of intent. When
the absolute_altitude_error333absolute_altitude_error is a derived fluent
becomes less than 1000 ft, the Altitude Alerting System shall initiate the
capture_alert and altitude_alert within 1 second.

The set of constraints, in s(CASP), in the body of

    happens(init_cap_alert)} capture the notion that the event has to occur within one second of the error dropping below 1000 ft. This shows the ease of expressing continuous real-valued constraints in s(CASP).
    % \textit

When the altitude error drops below 1000 ft, the alerting system shall raise
the capture alert within 1 second.

⬇

happens(init_cap_alert,T) :-

error_becomes_lt(1000,T1),T1#<T,T#<T1+1.

happens(init_alt_alert,T):-

error_becomes_lt(1000,T1),T1#<T,T#<T1+1.

initiates(init_cap_alert,cap_alert_on,T).

terminates(init_cap_alert,cap_alert_off,T).

initiates(init_alt_alert,alt_alert_on,T).

terminates(init_alt_alert,alt_alert_off,T).

REQ2: When the absolute_altitude_error becomes greater than 200 ft, the
Altitude Alerting System shall initiate the departure_alert and altitude_alert
within 1 second.

⬇

happens(init_dep_alert,T) :-

error_becomes_gt(200,T1),T1#<T,T#<T1+1.

initiates(init_dep_alert,dep_alert_on,T).

terminates(init_dep_alert,dep_alert_off,T).

happens(init_alt_alert,T) :-

error_becomes_gt(200,T1),T1#<T,T#<T1+1.

REQ3: When the absolute_altitude_error becomes either less than 180 ft or
greater than 1020 ft, the Altitude Alerting System shall terminate all
altitude_alerts within 1 second.

⬇

happens(term_all_alerts, T) :-

error_outside(180,1020,T1),T1#<T,T #<T1+1.

initiates(term_alerts,dep_alert_on,T).

terminates(term_alerts,dep_alert_off,T).

initiates(term_alerts,cap_alert_on,T).

terminates(term_alerts,cap_alert_off,T).

initiates(term_alerts,alt_alert_on,T).

terminates(term_alerts,alt_alert_off,T).

REQ4: While the altitude_selection_knob has not been adjusting within the
previous 5 seconds, the Altitude Alerting System shall determine
absolute_altitude_error as the absolute difference between the
selected_altitude and the barometric_altitude, otherwise
absolute_altitude_error shall be 0.

⬇

holdsAt(altitude_abs_error(E),T) :-

holdsAt(barometric_alt(B),T), holdsAt(selected_alt(S),T),

abs(S,B,E),T1 #= T-5, not happensIn(adjust(Value),T1,T).

holdsAt(altitude_abs_error(0), T) :-

happens(adjust(Value), T1, T1), T #>= T1, T #=< T1 + 5

happensIn(Event, T1, T2) :- happens(Event, T), T1 #=< T, T #=< T2.

REQ5: Upon initialization, the Altitude Alerting System shall consider the
altitude_selection_knob to be adjusting.

⬇

initiallyP(altitude_adjusting).

initiates(reset, altitude_adjusting, T).

###  4.2 Simulation of Alerting System Runs from Initial Conditions

The coding of the system requirements in event calculus and ASP is a logic
program that can be executed on the s(CASP) system. We can run various queries
to check if certain properties (e.g., safety properties) that we expect to
hold, are entailed by the encoding. We can construct various scenarios
(narratives, in event calculus parlance) and check that our requirements
specifications are indeed consistent, i.e., when we do simulation runs, we do
obtain a model. For example, we assume the following initial state of the
system and compute when the alerts will be raised (initiallyP(F) means that
initially fluent F is true and initiallyN(F) means initially fluent F is
false). The fluent start is used to define the rate of change of barometric
altitude as an event calculus trajectory444We have given encoding of
trajectory in Basic EC, the encoding of trajectory changes only slightly in
Extended EC. With appropriate narratives555We excluded `f`or lack of space as
part of the narrative given for the altitude alerting system, its end-to-end
behavior can be understood. To find out the time at which departure_alert_on
fluent holds, we issue the query Q1:

    ?- holdsAt(dep_alert_on, T).}
    which yields the answer \texttt

T > 30. This is because the altitude error becomes greater than 200 ft at T =
30.  ⬇ initiallyP(cap_alert_off).  initiallyP(dep_alert_off).
initiallyP(alt_alert_off).  initiallyP(altitude_adjusting).  ⬇
initiallyP(selected_alt(32000)).  initiallyN(selected_alt(V)) :- V \= 32000\.
trajectory(start,T1,barometric_alt(B),T2) :-  B #= 32200 + (T2 - T1) * 10\.

###  4.3 Adequacy of Requirements using Knowledge of Environmental Upsets

We can also test the behavior of the above scenario augmented with domain
knowledge. Our goal is test the adequacy of the requirements using prior
knowledge of environmental upsets. Adequacy means that in the presence of
known failures, the requirements are still consistent, i.e., a model exists
for them. In this case, the knowledge is simple and relates to Single Event
Upsets (or SEUs) [19] which can induce random resets of avionic software
systems. We can model the destructive power of the SEU using a reset event.
The reset event overrides the constructed internal state of the system,
forcing the initialization state to be re-established. Our goal in introducing
this reset is to ‘move the designer to consider how robust the system
initialization logic is to such transient resets. In ASP, an extraneous event
(e.g., reset) is represented as an abducible (a reset may or may not happen).
We would want to know if our system will still behave correctly in presence of
this extraneous event. The knowledge assumes that the reset event can only
override the cyber system (software state) and does not affect the continuous
state of the physical world. In case of the altitude alerting system, the
reset event can only affect the alert status and the altitude adjusting
status, but cannot affect the barometric altitude of the aircraft or the
selected altitude value. This in turn maps to the reset event terminating only
the fluents that describe the internal software state of the system. The reset
will not terminate fluents that describe the state evolution of continuous
quantities involving the real world. In general, fluents that are outside the
scope software system boundary are not affected by reset. The following
generic rule specifies that the reset event terminates every fluent associated
with the cyber system and initiates it to an initial value.

⬇

terminates(reset, Fluent, T) :- internal_state(Fluent).

initiates(reset, Default, T) :- default(Fluent, Default),
internal_state(Fluent).

Illustrating a Failing Run: Let the predicate $in\\_alerting\\_zone(T)$ denote
that the altitude error conditions for altitude alert ($REQ2$) are satisfied.
We show that the system fails to issue critical alerts in the presence of an
SEU-induced reset. We can check whether the altitude alerts are missed while
in an alerting zone using the query

    ?- in_alerting_zone(T), holdsAt(altitude_alert_off, T)} . With correct assumptions, the query would return ‘no models’ in s(CASP).
    By assuming an incorrect requirement in place of $REQ5$, the above query returns a binding in s(CASP), thereby showing a missed alert.   The incorrect requirement, $REQ5’$, is stated as: \textit

Upon initialization, the Altitude Alerting System shall consider the
altitude_selection_knob to be not Adjusting. More than likely, $REQ5$ may not
even be stated by a system designer, and if stated, it is likely that it will
be in the form of $REQ5^{\prime}$, as one would assume that the pilot would
not worry about adjusting the altitude knob during plane’s takeoff. Assume
that the plane is cruising at 32,200 feet at time $t$ and selected altitude
value is at 32300 ft. The plane’s barometric altitude increases at the rate of
10 ft per second. At $t+30$, the altitude error becomes greater than 200 ft.
Normally, the system should issue altitude alert at $t+30$. However, a reset
occurs at $t+30$. The reset sets the altitude to not adjusting as per
$REQ5^{\prime}$ and also turns the altitude alert off. Due to $REQ5^{\prime}$,
the altitude is treated as not adjusting for any time $T>t+30$. Thus, although
$in\\_alerting\\_zone(T)$ is true for $T>t+30$, the departure alert is never
turned back on. Therefore, an alert that should be issued to the pilot is
missed by the system. On the other hand, when $REQ5$ is used, the altitude
error is forced to be zero upon reset. This ensures that at $T=t+35$, (5
seconds from $t+30$) the error becomes 200 ft, thereby making
$in\\_alerting\\_zone(t+35)$ true. Again, per $REQ2$, the altitude alert would
be issued and no alert is missed. Adequacy of requirements with reset is
checked easily using query (Q2):

    ?-happens(reset,T),T1#>=T+5,in_alerting_zone(T1),holdsAt(alt_alert_off,T2),T2#>T1.
    }
    %\vspace

-0.3in The above query (Q2) asks s(CASP), that if in the event of a reset, is the altitude alert turned on after 5 seconds? If Q2 succeeds, then the requirements are adequate. If Q2 fails (say when we use $REQ5$ instead of $REQ5^{\prime}$) then the requirements are not adequate. This example shows the ease of validating requirements using EC and s(CASP) while using knowledge of SEUs. In summary, our ASP-coded EC model of the requirements are used to check for consistency and that they satisfy certain properties much in the spirit of model checking. This ASP and EC model can also be elegantly used for checking adequacy of requirements against known failure patterns. The failure patterns can be induced into the model through the use of abduction directly supported in s(CASP). This is possible due to near-zero “semantic gap” between the requirements and their mapping to ASP and EC. Checking adequacy is not direct in other approaches based on automata and Kripke structures.

###  4.4 Performance

The performance of s(CASP) for above queries is given in Table 1. The queries
were run on a quad code Intel i7-10510U processor at 1.8Ghz and 8GB RAM. For
comparison, we also attempted to run the queries on CLINGO [9] that is based
on grounding and SAT-solving. Due to large size of grounded program produced
from discretization of time and altitude values, the queries generally out of
memory. Note that Q3 is Q2 but rerun with REQ5 replaced by REQ5’.

Query | s(CASP) | CLINGO
---|---|---
Query Q1, Normal Run | 0.930s | ¿ 40 min
Query Q2, With Reset | 2.077s | ¿ 40 min
Query Q3 With Reset | 0.081s | ¿ 40 min
Table 1: Query Execution Times

##  5 Conclusions and Future Work

This paper presents a novel, systematic approach to formalizing requirement
modeling, which enables the integration and leveraging of domain knowledge
toward a computer-aided system validation, verification and assurance. We see
this work as a first step towards a new generation of Knowledge-Assisted
System Engineering (KASE) tooling. Our goal with this class of tooling is to
significantly reduce system development costs by systematically addressing
requirements specification defects much earlier in the development process.
Our approach is unique. It enables a formal examination of preliminary
conceptual models of intended behaviours, to assess adequacy and consistency
with respect to world assumptions, anticipated failure modes, and
environmental conditions. Our approach further enables the application of
domain knowledge and prior lessons learned (known deficiency areas) to be
fused with emerging specifications to more rapidly mature design robustness.
The primary catalyst for our success is not only the simplicity and elegance
of the EC, ASP, and the s(CASP) system, which allows us to approach human
levels of reasoning, it is also the intuitive formalism of EC being hidden
behind the CLEAR natural language style of specification, that enables the
power of formal methods to be introduced with minimal training and required
specialist knowledge. We hope this will remove some of the barriers to broader
spread of formal method adoption [7]. To the best of our knowledge, the work
reported here is novel, as the event calculus has not been used outside of AI.
Its application to system assurance here is novel. Even within AI, EC has not
been widely adopted because of the lack of sound, query-driven implementations
of negation-as-failure that can handle time faithfully as a continuous
quantity. Our s(CASP) implementation of ASP solves both these problems. The
work reported here also has practical value, as its results have benefited
engineering teams at a major aerospace company. Compared to model checking
techniques, our models are closer to the requirement specification, as states
do not have to be discerned in advance. Moreover, reasoning over continuous
time in EC is realized directly in the s(CASP) system. There are many
potential avenues for follow-on areas of research. First, we have already
internally established the feasibility to automatically generate the ASP
formal model from the MIDAS model and CLEAR requirements, using various
methods, including simple definite clause grammars to more state-of-the-art
symbolic natural language processing approaches [2]. We are also exploring
formal architectural reasoning using models of defined functional intent and
architectural intent with the abductive inference of component failures, as
well as modeling more complex cyber physical systems using the method reported
in this paper. Finally, we plan to apply our techniques to requirements
specification analysis of more systems, in a manner similar to UPPAAL
(https://uppaal.org/casestudies/).

## References

  * [1]
  * [2] Kinjal Basu et. al (2021): _Knowledge-driven Natural Language Understanding of English Text and its Applications_.  In: Thirty-Fifth AAAI Conference on Artificial Intelligence, AAAI 2021, 2021, AAAI Press, pp. 12554–12563.
  * [3] Joaquín Arias, Manuel Carro, Elmer Salazar, Kyle Marple & Gopal Gupta (2018): _Constraint answer set programming without grounding_.  TPLP 18(3-4):337-354, 10.1017/S1471068418000285.
  * [4] Joaquín et al. Arias (2019): _Modeling and Reasoning in Event Calculus Using Goal-Directed Constraint Answer Set Programming_.  In: LOPSTR’19, Springer, pp. 139–155, 10.1007/978-3-030-45260-5_9.
  * [5] C. Baral (2003): _Knowledge representation, reasoning and declarative problem solving_.  Cambridge University Press, 10.1017/CBO9780511543357.
  * [6] Devesh Bhatt, Brendan Hall et al. (2018): _The CLEAR Way To Transparent Formal Methods_.  In: 4th Workshop on Formal Integrated Development Environment (F-IDE), FLoC 2018.
  * [7] Jennifer Davis et al. (2013): _Study on the barriers to the industrial adoption of formal methods_.  In: International Workshop on Formal Methods for Industrial Critical Systems, Springer, pp. 63–77, 10.1007/978-3-642-41010-9_5.
  * [8] Dov Dori (2011): _Object-process methodology_.  In: Encyclopedia of Knowledge Management, Second Edition, IGI Global, pp. 1208–1220, 10.4018/978-1-59904-931-1.ch116.
  * [9] Martin Gebser et al. (2011): _Potassco: The Potsdam answer set solving collection_.  Ai Communications 24(2), pp. 107–124, 10.3233/AIC-2011-0491.
  * [10] M. Gelfond & Y. Kahl (2014): _Knowledge representation, reasoning, & design of intelligent agents: The answer-set programming approach_.  Cambridge Univ. Press, 10.1017/CBO9781139342124.
  * [11] B. Hall, D. Bhatt et al. (2018): _A CLEAR Adoption of EARS_.  In: IEEE EARS Workshop, pp. 14–15, 10.1109/EARS.2018.00010.
  * [12] B. Hall, J. Fiedor & Y Jeppu (2020): _Model Integrated Decomposition and Assisted Specification (MIDAS)_.  In: INCOSE Int’l Symp., 30(1), Wiley, pp. 821–841.
  * [13] G. H. Harman (1965): _The Inference to the Best Explanation_.  The Philosophical Review 74(1), pp. 88–95, 10.2307/2183532.
  * [14] Alistair Mavin & Philip Wilkinson (2010): _Big ears (the return of” easy approach to requirements engineering”)_.  In: 2010 18th IEEE International Requirements Engineering Conference, IEEE, pp. 277–282, 10.1109/RE.2010.39.
  * [15] Alistair Mavin et al. (2009): _Easy approach to requirements syntax (EARS)_.  In: 2009 17th IEEE International Requirements Engineering Conference, IEEE, pp. 317–322, 10.1109/RE.2009.9.
  * [16] Alistair Mavin et al. (2016): _Listens learned (8 lessons learned applying EARS)_.  In: 2016 IEEE 24th International Requirements Engineering Conference (RE), IEEE, pp. 276–282, 10.1109/RE.2016.38.
  * [17] Patrice Micouin (2014): _Model Based Systems Engineering: Fundamentals and Methods_.  John Wiley & Sons, 10.1002/9781118579435.
  * [18] Erik T. Mueller (2014): _Common Sense Reasoning: An Event Calculus Based Approach (2nd Edition)_.  Morgan Kaufmann.
  * [19] Eugene Normand (1996): _Single-event effects in avionics_.  IEEE Transactions on nuclear science 43(2), pp. 461–474, 10.1109/23.490893.
  * [20] M Sergot & R Kowalski (1986): _A logic-based calculus of events_.  New Generation Computing 4(1), pp. 67–95, 10.1007/BF03037383.
  * [21] Murray Shanahan (1999): _The event calculus explained_.  In: Artificial intelligence today, Springer, pp. 409–430, 10.1007/3-540-48317-9_17.
  * [22] WikePedia: _2015 Seville Airbus A400M crash_.  https://en.wikipedia.org/wiki/2015_Seville_Airbus_A400M_crash.