11institutetext: University of California Berkeley, Berkeley, CA 94720, USA
11email<EMAIL_ADDRESS>

# Explainable Automatic Grading with Neural Additive Models

Aubrey Condor 11 0000-0002-9482-7003 Zachary Pardos 11 0000-0002-6016-7051

###### Abstract

The use of automatic short answer grading (ASAG) models may help alleviate the
time burden of grading while encouraging educators to frequently incorporate
open-ended items in their curriculum. However, current state-of-the-art ASAG
models are large neural networks (NN) often described as "black box",
providing no explanation for which characteristics of an input are important
for the produced output. This inexplicable nature can be frustrating to
teachers and students when trying to interpret, or learn from an
automatically-generated grade. To create a powerful yet intelligible ASAG
model, we experiment with a type of model called a Neural Additive Model that
combines the performance of a NN with the explainability of an additive model.
We use a Knowledge Integration (KI) framework from the learning sciences to
guide feature engineering to create inputs that reflect whether a student
includes certain ideas in their response. We hypothesize that indicating the
inclusion (or exclusion) of predefined ideas as features will be sufficient
for the NAM to have good predictive power and interpretability, as this may
guide a human scorer using a KI rubric. We compare the performance of the NAM
with another explainable model, logistic regression, using the same features,
and to a non-explainable neural model, DeBERTa, that does not require feature
engineering.

###### Keywords:

Explainable AI Automatic Grading ASAG Neural Additive Models.

##  1 Introduction

It has been shown that the use of open-ended (OE) items is beneficial for
student learning due to the generation effect [3] or in combination with self-
explanation [6]. However, assessing OE items is time consuming for teachers
[16] and consequently, educators default to using multiple choice (MC)
questions instead. Automatic short answer grading (ASAG) may alleviate this
time burden while encouraging educators to frequently incorporate OE items in
their curriculum.

Many of the best performing ASAG models include some variation of a deep
neural network (NN) [15]. NNs are impressive predictors for high dimensional
inputs like text embeddings, but predictive power tends to come at the cost of
intelligibility. These models are often described as "black box" which means
that users only have access to inputs and outputs, yet no information as to
the process in-between. Further, unlike explainable models such as a logistic
regression (LR), with NN ASAG models, we are not choosing which features the
model may consider, and the model provides no explanation for which
characteristics of an input are important for the produced output. Whereas
educators typically use a scoring rubric to provide an explanation to
stakeholders about which features of a given response substantiate the
assigned score, and despite that attempts have been made to integrate scoring
rubrics into ASAG models [8], NN ASAG models provide no justification for
their predictions. The inexplicable nature of ASAG models can be frustrating
to both teachers and students when trying to make sense of, or learn from an
automated grade. Teachers are unable to monitor student understandings at a
fine-grained level, and students may not productively learn without knowing
why they received a certain grade. It is critical that researchers find more
explainable models if ASAG systems are to be of greater practical use.

In an effort to create a powerful yet intelligible ASAG model, we experiment
with a type of model that combines the performance of a NN with the
explainability of an additive model called a Neural Additive Model (NAM) [1].
NAMs allow us to visually examine the contribution of each feature to the
final predicted score for each response, similar to testing the significance
of a regression coefficient. For this type of model, we must engineer features
of a response as inputs to the model instead of allowing the model to create
its own text features like typical Large Language Model (LLM) classifiers do.
The research questions we seek to answer include, (1) can NAMs provide
intelligible automatic grading such that stakeholders can understand which
features of a response are important for its prediction, and (2) is the
predictive performance of NAMs better than that of legacy explainable models
like a LR and commensurate with that of an LLM classifier? This research is
unique as NAMs have not yet been explored for educational applications, or
more specifically for explainable automatic grading.

We demonstrate our NAM approach for ASAG with one item designed under a
Knowledge Integration (KI) perspective from the learning sciences and
corresponding rubric to guide our feature engineering. KI is a framework for
strengthening science understanding that emphasizes incorporating new ideas
and sorting out alternative perspectives with evidence [22]. We hypothesize
that the inclusion (or exclusion) of predefined KI ideas as features will be
sufficient for the NAM to have good predictive power, as this is precisely
what would guide a human scorer using the KI rubric. We compare the
performance of the NAM with another explainable model, a LR, using the same
features as those used for the NAM, and to a non-explainable model, DeBERTa,
an improved version of the popular BERT model [17]. The DeBERTa model does not
utilize the engineered features that are used by the NAM or LR model, but
creates its own features from the text. For a more extensive comparison of the
performance of these three models, we extend our analysis to include results
from five additional questions from an open-source data set that has been used
extensively in ASAG research - the Automatic Student Assessment Prize (ASAP)
Short Answer Scoring (SAS) data [28].

##  2 Related Work

In this section, we describe previous works relating to Explainable AI and
more particularly Explainable ASAG, as well as applications of NAMs.

###  2.1 Explainable AI and ASAG

A recent review article outlines the increase in demand for explainable AI
(XAI), especially for people who are affected by AI driven decisions [36].
They describe the rising popularity of NN models for their state-of-the-art
performance, despite that their inference processes are not known or
interpretable. Many researchers aim to increase the explainability of AI
systems. Initiatives such as DARPA’s "Explainable AI (XAI) program" [14] and
the European Union’s "General Data Protection Regulation" which demands
citizen’s rights to an explanation for decisions made by AI [13] are
contributing to an increased demand for XAI.

Only recently are researchers beginning to think about XAI in terms of ASAG
models. Schlippe et al. [31] surveyed over 70 educators about preferences for
explainability in ASAG and learned that most prefer to see matches between
student answers and exemplary answers to justify scores over other
explainability methods like highlighting the importance of certain words.
Poulton et al. [27] proposed an XAI framework using SHAP values to assess
popular LLMs for ASAG such as BERT and RoBERTa. Zeng et al. [37] investigated
whether autograding models align with human graders in terms of the important
words they use when assigning a grade by conducting a randomized controlled
trial to see if highlighting words deemed important by an autograder can
assist human grading. Singh et al. [32] introduced Summarize and Score (SASC)
to explain ‘text modules’, which map text to scalar values within LLMs,
providing a natural language explanation of the module’s selectivity. Finally,
Tornqvist et al. [33] introduced ExASAG, an explainable framework for ASAG
that generates natural language explanations for predictions, and Condor and
Pardos train a reinforcement learning agent to alter students’ short responses
so that these alterations provide insight as to why a ASAG score was assigned
[7].

###  2.2 Applications of Neural Additive Models

Although Neural Additive Models (NAMs) were introduced only a few years ago
[1], there has been substantial interest in their use as the machine learning
community has increased efforts to promote understandable models. Some
researchers have proposed altered or improved versions of NAMs for specific
applications. Mariotti et al. [24] explored the tension between
interpretability and performance of NAMs and introduce a constrainable NAM
(CNAM) which includes specifications for model regularization. The CNAM model
is able to consider the performance-interpretability tradeoff during training,
and is demonstrated to work well for both regression and classification tasks
[24]. Luber et al. [23] introduced Structural NAMs (SNAMs) which incorporate
neural splines to a typical NAM. SNAMs offer improved intelligibility over
NAMs by enabling direct interpretation of estimated parameters, as well as
quantification of parameter uncertainty and post-hoc analysisr. Bouchiat et
al. [4] take a Bayesian perspective and develop a Laplace approximated NAM
(LA-NAM) to enhance interpretability, and Jo and Kim [18] introduced NAMs for
multivariate time series modeling called NAM nowcasting (NAM-NC). Others have
used NAMs for applications in finance [5], mortality prediction [26], and
survival analysis [34]. Importantly, we found no previous works that used NAMs
for educational applications.

##  3 Background

In this section, we provide a description of the data, and a brief overview of
NAMs, LR, and the DeBERTa model.

###  3.1 The Data

The KI item used for this project was collected during a previous research
project at the Web-based Inquiry Science Environment (WISE) research center at
the University of California, Berkeley consisting of OE science items designed
for middle school students [30]. Students accessed the items via an online
classroom system, and responses were scored with a Knowledge Integration (KI)
rubric. For this project, we use one item from a unit about the physics of
sound waves that engages students to refine their ideas about concepts such as
wavelengths, frequency and pitch. For this item, students must distinguish how
the pitch of sound made by tapping a full glass of water compares to the pitch
made by tapping an empty glass. They are asked to explain why they think the
pitch of the sound waves may be the same or different for the two glasses
[30]. See Figure 1 for a visual of the Soundwaves item. The KI data include
1,313 OE student responses that were carefully rated from 1 to 5 by subject
matter experts. A detailed rubric was created by researchers and used for
scoring which includes a description of the rating level, examples of
correct/incorrect mechanisms and conclusions, and exemplar student responses
that would fall into each category. 2 provides a sample of the KI scoring
rubric.

The ASAP data used to extend our model comparison is from a 2012 Kaggle
competition sponsored by the Hewlett Foundation, and consists of almost 13,000
short answer responses to 10 science and English questions [28]. We used only
the 5 science questions to match the domain of the KI data. Each of the five
items have between 1300-1800 OE student responses and corresponding human-
rated scores between 1 and 4 with 4 being the most correct score. Each item
has an expert created scoring rubric, and we aggregate the 5 items for
training, testing and the resulting model comparisons.

Figure 1: The Sound Waves item bundle Figure 2: The KI Scoring Rubric

###  3.2 Neural Additive Models

Neural Additive Models (NAMs) impose a restriction on the structure of NNs in
order to make the model interpretable. NAMs belong to a family of models
called Generalized Additive Models (GAMs) which have the form:

$g(E[y])=+f_{1}(x_{1})+f_{2}(x_{2})+...+f_{k}(x_{k})$

Where $x=(x_{1},x_{2},\dots,x_{k})$ represents the input with K features and
each $f_{i}$ is a univariate shape function with $E[f_{i}]=0$ [1]. NAMs learn
a linear combination of jointly-trained NNs where each NN attends to a single
input feature. They can approximate complex, high-dimensional functions and
their predictions are easily interpretable. NAMs provide advantages over
standard GAMs such as their superior scalability with the use of GPU/TPU
hardware developed for NNs, differentiability, and visual interpretability
[1].

###  3.3 Logistic Regression

Multinomial LR is a classification model that predicts probabilities of
different outcomes for a categorical dependent variable. To generalize to a
K-class setting, the model runs K-1 independent binary LR models where one
outcome is chosen as a “pivot” and other K-1 outcomes are separately regressed
against the pivot outcome. We use the Limited-memory Broyden-Fletcher-
Goldfarb-Shannon (LBfGS) algorithm for optimization [12], and incorporate L2
regularization.

###  3.4 DeBERTa

DeBERTa (Decoding-enhanced BERT with disentangled attention) is a Large
Language Model (LLM) that improves upon the popular BERT [10] model by using a
different version of the standard attention mechanism [35] and a novel
position encoding method. DeBERTa uses a ‘disentangled attention’ mechanism
where each word is represented with two vectors - one to encode the word’s
context, and another to encode the word’s relative position. Further, an
enhanced masked decoder incorporates absolute positions to predict masked
tokens during training. Finally, DeBERTa uses a unique adversarial training
method to fine-tune the model’s generalizability performance on downstream
tasks. We use an improved yet smaller version of the original DeBERTa model
called DeBERTaV3-base [17], consisting of 12 layers, a hidden size of 768, and
86 million parameters.

##  4 Methods

To fit a NAM, we utilize the python library, "nam" [19], created by the
researchers who introduced NAMs for easy implementation of the model. We
altered the package code to accommodate our multi-class classification setting
and allow for our evaluation procedure and corresponding metrics. Further
described in the "Evaluation" section below, we use a cross validation (CV)
evaluation procedure for model comparison, so for each model class, we train a
total of 10 models on a different sample of the data. We train each of the 10
NAMs for 120 epochs with a batch size of 64 and a learning rate of 0.002. An
epoch is one complete pass of the training data through the model, batch size
represents the number of samples used in a single forward and backward pass
through the network, and the learning rate governs the pace at which the model
updates its parameter estimates. Further, we use a regularization technique
called dropout with a value of 0.15 to help avoid over fitting.
Hyperparameters were chosen using a standard grid search where different
combinations of values are tested on a held-out validation set, and those that
result in the best performance are chosen. Each of the DeBERTa models were
trained for six epochs using a batch size of eight, and a learning rate of
0.0007. The LR models were fit using the scikit-learn "LogisticRegression"
python library.

###  4.1 Feature Engineering

To create features from student response texts for the NAM and LR models, we
utilize KI phrases from the rubric that represent both correct, and incorrect
mechanisms and conclusions. We chose to use KI ideas to create features
because these are exactly what a human identifies to assign a grade to a
student response within the KI framework. Because students often express their
ideas using a variety of words (i.e., the idea that "sound moves faster in
air" could be stated as "sound travels faster through air"), we scan the
responses for n-grams that are similar to the chosen KI phrases. An n-gram is
a series of n consecutive words within a string of text. For example, 2-grams
of the phrase "the dog ran" would be "the dog" and "dog ran". We use n-grams
of size n=1 through 5 for each response as students can also express the same
ideas in a different number of words. Features for the ASAP science questions
were created in the same manner (and use the same number of features for each
item), although we provide the feature phrases only from the KI data for
demonstration.

To calculate the similarity between a given feature phrase and an n-gram from
a student response, we first embed both the n-gram and the phrase using
sentence-BERT. Sentence-BERT utilizes siamese and triplet network structures
to create semantically meaningful sentence embeddings [29]. We then calculate
the cosine similarity of the KI phrase and the n-gram. Cosine similarity is a
metric used to measure the similarity of text using the cosine of the angle
between the vector embeddings of the two texts being compared. Once the KI
phrase has been compared with all respective n-grams in a response, we use the
largest cosine similarity found as the feature for that phrase. Thus the
feature representation for each phrase provides a measure of whether or not
some form of the phrase is included in the student response. We use 34 phrases
and 28 key words from the KI Sound Wave Rubric (shown in Figure 2) for a total
of 62 features. The same number of phrases and words were chosen for each ASAP
item from their corresponding scoring rubrics. The KI phrases are shown below,
and the individual words were chosen from the phrases:

"I don’t know", "idk", "sound bounces", "water blocks sound", "water mutes
ringing", "sound moves more in air", "sound moves less in water", "sound
echoes", "sound sinks in water", "sound moves fast in air", "frequency is
height of wave", "pitch higher in full glass", "the density is different",
"water is more dense", "water vibrates less", "affects vibration", "pitch is
lower in water", "pitch is different", "higher frequency in air", "frequency
is different", "sound moves faster in water", "water has more mass", "mass is
different", "sound is denser in water", "sound is slower in water", "amplitude
is number of waves", "pitch lower in empty glass", "air is less dense", "empty
glass vibrates more", "vibration is different", "pitch is higher in air",
"lower frequency in water", "sound moves slower in air", "empty glass less
mass"

###  4.2 Evaluation

We use a standard ASAG evaluation metric in our results table: Quadratic
Weighted Cohen’s Kappa (QWK). QWK reports the agreeability between two scores
beyond random chance - a more robust measure than accuracy. A 5x2 cross
validation (CV) paired t-test was used to evaluate the statistical
significance of the difference in models. The 5x2 CV paired t-test is based on
five iterations of twofold cross-validation, and is presented in [11] as the
recommended approximate statistical test for whether one machine learning
model outperforms another because of it’s more acceptable type I error, and
stronger statistical power than other methods such as McNemar’s test, or a
paired t-test based on 10-fold CV. We used the QWK metric for the 5x2 CV
paired t-tests (with 5 degrees of freedom). We perform a t-test comparing the
QWK metric for the NAM to the LR model and the DeBERTa model. The null
hypothesis for each t-test states that the QWK metric for the NAM is no
different than the QWK metric for LR and DeBERTa. Further, we report the
magnitude of the difference in performance of the NAM to LR and DeBERTa using
the QWK metric from the CV model fits with a Cohen’s D effect size. Cohen’s D
is essentially a difference of means, scaled by a standard deviation value
[20]. A low value of Cohen’s D which represents a small magnitude of
difference is typically below 0.20, and a large value would be around 0.80 or
higher. We use a pooled standard deviation to calculate Cohen’s D, which
represents a weighted combination of the standard deviations from each model.

##  5 Results

Comparisons of each ASAG model’s performance are presented in tables 1 and 2.
In Table 1, we present results of the 5x2 CV paired t-test using the QWK
metric for both the KI data and the aggregate ASAP data. We compare the NAM
against both LR and DeBERTa. For the KI data, the NAM performs better than LR
at a significance level of 0.05 (t=2.1732, p=0.0409), and the NAM performs
slightly worse than the DeBERTa model but not statistically at a significance
level of 0.05 (t=-0.8410, p=0.2194). Further, for the ASAP data, the NAM does
not perform significantly better than the LR at a significance level of 0.05
(t=1.691, p=0.0758) although the p-value is quite small. For the ASAP data, we
also see that the NAM performs worse than the DeBERTa model but not at a
statistically significant level of 0.05 (t=-1.741, p=0.071), although
similarly, the p-value is quite small. As shown by the CV average metrics in
Table 2, on average, DeBERTa performs better than the NAM and the NAM performs
better than the LR across both data sets. Additionally, in Table 1 we present
the Cohen’s D effect sizes for the difference in performance from the NAM to
the two other models. The magnitude of difference between the NAM and both
models is quite large for the KI data, with effect sizes greater than 1.0 for
each comparison. For the ASAP data, the effect size is moderate when comparing
the NAM and LR.

Table 1: 5x2 Cross Validation Paired t-test (5 df) and Cohen’s D Effect Sizes NAM versus: | t(5) | p-val | Cohen’s D
---|---|---|---
| KI data |  |
DeBERTa | -0.8410 | 0.2194 | -1.111
Log Reg | 2.1732 | 0.0409 | 1.003
| ASAP data |  |
DeBERTa | -1.741 | 0.0711 | -1.070
Log Reg | 1.691 | 0.0758 | 0.222

A visual of feature importance from the NAM for the KI data is presented in
Figure 3. We show the top 40 most important features in order of increasing
importance. We see that the feature importance for the words "density" and
"higher" is much larger than that of most other words and phrases. The phrase
"the density is different" also has a comparatively high importance. Further,
in Figure 4, we present NAM Shape functions for the top 8 important words and
phrases. As each feature of the NAM is handled independently by a learned
shape function, we can see how the model makes its predictions, or in our case
ratings, by graphing the shape functions for each feature. The scores are
averaged and centered so that the shape functions can be directly compared on
the same scale. As our NAM is a multi-class classifier, we can visualize five
shape functions for each feature representing each of the five classes. We
choose to show only the shape function for the lowest rating class and the
highest rating class to retain the most important information yet avoid a
messy, undecipherable visual. However, Figure 5 provides an example of shape
plots with all five shape functions shown for reference. A negative y-axis
score signals a low probability of a certain class, and a positive score
represents a high probability. Noticeably, the shape functions are jagged-like
with sharp jumps. [1] emphasizes the unique capability of NAMs to model highly
“jumpy” 1-dimensional functions due to the use of an exponential-centered
(ExU) non-linear function, whereas standard neural network models are bias
towards smoothness with the use of Rectified Linear Unit (ReLU) functions. We
could hypothesize that the true shape functions for the features in this
project are more smooth, and in future work we could employ regularization
techniques such as weight decay to provide smoother shape functions. On the
same plot, we see the data density in the form of pink shaded bars. The darker
the shade of pink, the more data there is in that region. In some areas with a
really low data density, the model may not have had enough data to adequately
learn the shape function.

Table 2: 5x2 Cross Validation QWK Averages Model | KI data | ASAP data
---|---|---
DeBERTa | 0.7475 | 0.7176
Log Reg | 0.6929 | 0.6342
NAM | 0.7174 | 0.6473
Figure 3: NAM Mean Feature Importance

##  6 Discussion

The CV averages in Table 2 give us an idea for how well each model performs
for ASAG. We are not surprised to see that the NAM provides more predictive
power than the LR, but not as much as the DeBERTa model. The DeBERTa model
uses the entire student response as input text and learns many hundreds of
features for its predictions, whereas with the NAM and LR, we are limited to
the particular features we choose and it would be an unreasonable task to
engineer the same number of features as the DeBERTa model learns.
Additionally, the results of the 5x2 CV paired t-test in Table 1 give evidence
that the NAM outperforms the LR model - across both data sets at a
significance level of 0.10 and more notably with the KI data at a significance
level of 0.05. Importantly, the NAM and the LR model use the same exact
features. Further, the Cohen’s D effect size values show that the magnitude of
difference in performance of the models is quite large for the KI data, and
smaller but still notable for the ASAP data. Not only is the difference
statistically significant, but it seems to be practically meaningful.

Figure 4: NAM shape functions of the highest and lowest rating category for
the top 8 most important phrases/words. The y-axes show the log odds of
predicting a given rating category, and the x-axes represent the range of
similarity scores. The pink shades represent the data density at varying
similarity scores.

Figure 3 allows us to visually interpret the ASAG NAM through overall feature
importance, and Figure 4 gives us more detail about how each feature
contributes with NAM shape functions. Human grading of OE responses can often
consist of identifying which ideas - correct or incorrect - are included in a
student’s response. The KI scoring rubric allowed us to identify which ideas
(phrases and words) may help guide the NAM to correctly rate student
responses, and the resulting NAM visualizations enable us to discern how these
ideas contribute to an automated rating. The feature importance plot can help
stakeholders identify which ideas are most indicative of the student’s grade,
and the NAM shape function gives a notion of how the directionality of the
model’s prediction changes for different values of different rating classes.
For example, for the “higher” feature, the probability for the highest rating
class (shown by the green line in Figure 4) goes down significantly with
decreasing cosine similarities. So, including the word “higher” in a response
seems to correspond to higher scores. Additionally, with the density shading
in the NAM shape function plots, we can infer the density of responses that
include the given word or phrase. For example, the words “more”, and “higher”
both have dense regions near high cosine similarities, so we can assume that
many responses include these words or similar words. This can give an educator
insight into which ideas students are more or less likely to write about.

Figure 5: An Example of NAM shape functions with all rating categories

Limitations of using a NAM for ASAG include that the performance of the model
in terms of its match to human ratings is less than that of a LLM ASAG model.
Also, engineering features for the NAM is more costly than preparing student
responses for a LLM ASAG model, as the model will create features from the
text itself. Further, we only show results for science questions so we cannot
conclude that our results would generalize to other short answer items from
different subject areas. Our results suggest that the NAM may be a sufficient
alternative to legacy explainable models, like a LR. Not only did the NAM
generally perform better than the LR with the same features, but the NAM
provides easily interpretable visualizations of the model’s prediction
functions which can give an educator insights about their students’
understanding of the item content. Many researchers in the field of Learning
Analytics still use LR over NNs for classification when the advantages of
intelligibility outweigh that of performance [21][9][2][25]. The Learning
Analytics community may benefit from investigating the use of NAMs for various
classification tasks as well. In future work, we hope to experiment with using
NAMs for ASAG with different question types from different domains, and
perform qualitative interviews with educators to see if the NAM visualizations
are understandable and helpful.

## References

  * [1] Agarwal, R., Frosst, N., Zhang, X., Caruana, R., Hinton, G.E.: Neural additive models: Interpretable machine learning with neural nets. arXiv preprint arXiv:2004.13912 (2020)
  * [2] Alonso-Fernández, C., Martínez-Ortiz, I., Caballero, R., Freire, M., Fernández-Manjón, B.: Predicting students’ knowledge after playing a serious game based on learning analytics data: A case study. Journal of Computer Assisted Learning 36(3), 350–358 (2020)
  * [3] Bertsch, S., Pesta, B.J., Wiscott, R., McDaniel, M.A.: The generation effect: A meta-analytic review. Memory & cognition 35(2), 201–210 (2007)
  * [4] Bouchiat, K., Immer, A., Yèche, H., Rätsch, G., Fortuin, V.: Laplace-approximated neural additive models: Improving interpretability with bayesian inference. arXiv preprint arXiv:2305.16905 (2023)
  * [5] Chen, D., Ye, W.: Monotonic neural additive models: Pursuing regulated machine learning models for credit scoring. In: Proceedings of the Third ACM International Conference on AI in Finance. pp. 70–78 (2022)
  * [6] Chi, M.T., De Leeuw, N., Chiu, M.H., LaVancher, C.: Eliciting self-explanations improves understanding. Cognitive science 18(3), 439–477 (1994)
  * [7] Condor, A., Pardos, Z.: A deep reinforcement learning approach to automatic formative feedback. International Educational Data Mining Society (2022)
  * [8] Condor, A., Pardos, Z., Linn, M.: Representing scoring rubrics as graphs for automatic short answer grading. In: International Conference on Artificial Intelligence in Education. pp. 354–365. Springer (2022)
  * [9] Deho, O.B., Zhan, C., Li, J., Liu, J., Liu, L., Duy Le, T.: How do the existing fairness metrics and unfairness mitigation algorithms contribute to ethical learning analytics? British Journal of Educational Technology 53(4), 822–843 (2022)
  * [10] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [11] Dietterich, T.G.: Approximate statistical tests for comparing supervised classification learning algorithms. Neural Computation 10(7), 1895–1923 (1998). https://doi.org/10.1162/089976698300017197
  * [12] Fletcher, R.: Practical methods of optimization. John Wiley & Sons (2000)
  * [13] Goodman, B., Flaxman, S.: European union regulations on algorithmic decision making and a ”right to explanation”. AI magazine 38(2), 781–796 (2017)
  * [14] Gunning, D., Aha, D.: Darpa’s explainable artificial intelligence (xai) program. AI magazine 40(2), 44–58 (2019)
  * [15] Haller, S., Aldea, A., Seifert, C., Strisciuglio, N.: Survey on automated short answer grading with deep learning: from word embeddings to transformers. arXiv preprint arXiv:2204.03503 (2022)
  * [16] Hancock, C.L.: Implementing the assessment standards for school mathematics: Enhancing mathematics learning with open-ended questions. The Mathematics Teacher 88(6), 496–499 (1995)
  * [17] He, P., Gao, J., Chen, W.: Debertav3: Improving deberta using electra-style pre-training with gradient-disentangled embedding sharing. arXiv preprint arXiv:2111.09543 (2021)
  * [18] Jo, W., Kim, D.: Neural additive models for nowcasting. arXiv preprint arXiv:2205.10020 (2022)
  * [19] Kayid, A., Frosst, N., Hinton, G.E.: Neural additive models library (2020)
  * [20] Kelley, K., Preacher, K.J.: On effect size. Psychological methods 17(2), 137 (2012)
  * [21] Le, C.V., Pardos, Z.A., Meyer, S.D., Thorp, R.: Communication at scale in a mooc using predictive engagement analytics. In: Artificial Intelligence in Education: 19th International Conference, AIED 2018, London, UK, June 27–30, 2018, Proceedings, Part I 19. pp. 239–252. Springer (2018)
  * [22] Linn, M.C.: Designing the knowledge integration environment. International Journal of Science Education 22(8), 781–796 (2000)
  * [23] Luber, M., Thielmann, A., Säfken, B.: Structural neural additive models: Enhanced interpretable machine learning. arXiv preprint arXiv:2302.09275 (2023)
  * [24] Mariotti, E., Moral, J.M.A., Gatt, A.: Exploring the balance between interpretability and performance with carefully designed constrainable neural additive models. Information Fusion p. 101882 (2023)
  * [25] Misiejuk, K., Wasson, B., Egelandsdal, K.: Using learning analytics to understand student perceptions of peer feedback. Computers in human behavior 117, 106658 (2021)
  * [26] Moslehi, S., Mahjub, H., Farhadian, M., Soltanian, A.R., Mamani, M.: Interpretable generalized neural additive models for mortality prediction of covid-19 hospitalized patients in hamadan, iran. BMC Medical Research Methodology 22(1), 339 (2022)
  * [27] Poulton, A., Eliens, S.: Explaining transformer-based models for automatic short answer grading. In: Proceedings of the 5th International Conference on Digital Technology in Education. pp. 110–116 (2021)
  * [28] Prize, A.S.A.: The hewlett foundation: Automated essay scoring (2019)
  * [29] Reimers, N., Gurevych, I.: Sentence-bert: Sentence embeddings using siamese bert-networks. arXiv preprint arXiv:1908.10084 (2019)
  * [30] Riordan, B., Bichler, S., Bradford, A., King Chen, J., Wiley, K., Gerard, L., Linn, M.: An empirical investigation of neural methods for content scoring of science explanations. In: Proceedings of the fifteenth workshop on innovative use of NLP for building educational applications (2020)
  * [31] Schlippe, T., Stierstorfer, Q., Koppel, M.t., Libbrecht, P.: Explainability in automatic short answer grading. In: International Conference on Artificial Intelligence in Education Technology. pp. 69–87. Springer (2022)
  * [32] Singh, C., Hsu, A.R., Antonello, R., Jain, S., Huth, A.G., Yu, B., Gao, J.: Explaining black box text modules in natural language with language models. arXiv preprint arXiv:2305.09863 (2023)
  * [33] Tornqvist, M., Mahamud, M., Guzman, E.M., Farazouli, A.: Exasag: Explainable framework for automatic short answer grading. In: Proceedings of the 18th Workshop on Innovative Use of NLP for Building Educational Applications (BEA 2023). pp. 361–371 (2023)
  * [34] Utkin, L., Konstantinov, A.: An extension of the neural additive model for uncertainty explanation of machine learning survival models. In: Cyber-Physical Systems: Intelligent Models and Algorithms, pp. 3–13. Springer (2022)
  * [35] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, Ł., Polosukhin, I.: Attention is all you need. Advances in neural information processing systems 30 (2017)
  * [36] Xu, F., Uszkoreit, H., Du, Y., Fan, W., Zhao, D., Zhu, J.: Explainable ai: A brief survey on history, research areas, approaches and challenges. In: CCF international conference on natural language processing and Chinese computing. pp. 563–574. Springer (2019)
  * [37] Zeng, Z., Li, X., Gasevic, D., Chen, G.: Do deep neural nets display human-like attention in short answer scoring? In: Proceedings of the 2022 conference of the north American chapter of the Association for Computational Linguistics: Human language technologies. pp. 191–205 (2022)


# Modeling Uplink Coverage Performance in Hybrid Satellite-Terrestrial
Networks

Bassel Al Homssi, , and Akram Al-Hourani B. Al Homssi and A. Al-Hourani are
with the School of Engineering, RMIT University, Melbourne, Australia. E-mail:
<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
###### Abstract

Once deemed a far-fetched vision, emerging deployments of massive satellite
constellations will soon offer true global coverage. When these constellations
overlay the evolving terrestrial networks, a new hybrid continuum is formed
which has the potential to provide uninterrupted coverage. In this paper, we
provide an analytic framework for the uplink coverage probability in hybrid
satellite-terrestrial networks. The framework extends well-developed
terrestrial models that utilize tools from stochastic geometry by
incorporating an additional layer that fits the emerging next generation
satellite-terrestrial networks. The paper captures the impact of both (i) the
constellation size, and (ii) the terrestrial base station density on the
coverage of the uplink traffic which is dominant in applications relying on
wireless sensor networks, such as the Internet of Things. This framework
provides insights that guide the design of hybrid network infrastructures for
a desired quality of service.

###### Index Terms:

Stochastic geometry, IoT-over-satellite, massive satellite constellation,
uplink, LEO, hybrid network.

##  I Introduction

Next generation terrestrial networks are expected to provide improved
mobility, connectivity, and capacity to rival the current soar in demand for
higher quality and faster uplink communications, especially with the
deployment of massive wireless Internet-of-Things (IoT) sensor networks. A
large number of these IoT sensors are expected to be dispersed in remote
locations serving various applications such as smart agriculture, asset
tracking, mining, oil rigs, and environmental sensing. However, remote and
off-shore deployments are not typically supported by terrestrial coverage due
to technical difficulties and cost limitations. Therefore, it is imperative to
have a more flexible network architecture that incorporates both terrestrial
and satellite technologies to facilitate seamless coverage continuum. Driven
by this vision, the third generation partnership project (3GPP) is actively
looking into hybrid satellite-terrestrial networks under the 5G framework [1].
The satellite link can provide complementary coverage to the underlying
terrestrial infrastructure and ensure service continuity. Fig. 1 depicts the
concept of hybrid satellite-terrestrial network where a seamless coverage is
potentially attained [2].

Figure 1: An illustration of the hybrid satellite-terrestrial continuum, where
terrestrial gaps are bridged by overlaying satellite coverage.

The downside is that unlike terrestrial base stations (BS), satellites are
usually very expensive and require long lead time to deploy. This is only
exacerbated with the need for massive numbers of satellites required to
achieve adequate global coverage. As a result, network planning is of
paramount importance to sustain cost-effective deployment which meets operator
targets. Analytic approaches can usually provide deeper insights into the
different parameters impacting the performance of the complex satellite
constellation networks. Tools from stochastic geometry are widely deployed to
capture the performance of terrestrial networks, these tools have been
recently extended to analyze the behavior of satellite networks [3]. Example
of recent work can be found in [4, 5] for analyzing the coverage probability
by relying on the log-distance path-loss channel model typically used to
characterize terrestrial links. However, the vast proportion of the satellite
communication occurs under free-space conditions and only interacts with the
clutter layer near ground [3, 6]. Moreover, the literature lacks an analytic
framework that provides the uplink coverage for hybrid satellite-terrestrial
networks. Thus, in this letter we provide an analytic framework that models
the uplink coverage probability for hybrid satellite-terrestrial networks in
rural environments using satellite specific path-loss and shadowing model. The
hybrid coverage is assumed to be modeled using two independent and
homogeneously distributed sets of receivers; (i) a satellite constellation,
and (ii) terrestrial BS. Accordingly, the framework presents a relationship
between the coverage probability and the operation parameters for both
networks, including; (i) the number of satellites, (ii) base station (BS)
density, and (iii) ground devices density. The framework also presents an
expansion strategy based on a given target quality-of-service (QoS) and the
projected devices increase. The contributions of this work are summarized as
follows:

  * •

It presents an analytic approach for the uplink coverage probability in
satellite-terrestrial networks suitable for IoT over satellite applications.

  * •

It derives analytic relations between the number of satellites and BS density
for a prescribed QoS.

  * •

It depicts an expansion strategy for satellite-terrestrial networks in terms
of satellite numbers and BS density to cater for the projected increase in
ground users.

##  II System Model

Figure 2: A user-centric simplified representation of the hybrid network where
the user is jointly being served by a terrestrial BS and a satellite.

###  II-A Geometric Model

We define three independent point processes that respectively represent the
locations of (i) the satellite constellation, (ii) terrestrial BS, and (iii)
users. The satellite locations are considered to be randomly and homogeneously
scattered on a constellation sphere of altitude $h$ above the mean sea level,
each satellite follows a random circular orbit with the same radius. This
assumption preserves the uncorrelation and uniform distribution overtime. We
show later that this assumption provides the upper bound and behaves closely
in comparison to deterministic Walker satellite constellations. Therefore, the
satellite locations can be approximated by the homogeneous Poisson point
process (PPP) as shown in [3]. Consider a constellation is composed of
$N_{\mathrm{s}}$ satellites orbiting on a spherical surface with an altitude
of $h$, then the average satellite density is given as,

$\lambda_{\mathrm{s}}=\frac{N_{\text{s}}}{4\pi(R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}+h)^{2}},$
(1)

where $R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}$ is the Earth’s average
radius. Each user is assumed to connect to the satellite that is closest to
its local zenith due to better radio propagation conditions. Using geometric
reasoning, we represent the user-centered zenith angle $\varphi^{\prime}$ in
terms of the Earth-centered zenith angle $\varphi$, as follows,

$\varphi^{\prime}=\operatorname{acot}\left(\frac{\cos\varphi-\alpha}{\sin\varphi}\right),$
(2)

where
$\alpha=R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}/(R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}+h)$,
and both angles are illustrated in Fig. 2. Accordingly, the probability
density function (PDF) of the contact angle $\varphi_{\mathrm{o}}$ is obtained
in [3] as follows,

$f_{\varphi_{\mathrm{o}}}(\varphi)=\frac{N_{\mathrm{s}}}{2}\sin\varphi\exp\left(-\frac{N_{\mathrm{s}}}{2}[1-\cos\varphi]\right).$
(3)

On the other hand, the terrestrial BS locations are represented with an
independent and homogeneous PPP on Earth’s surface $\Phi_{\mathrm{b}}$ on
$\mathbb{S}^{2}$ with density ${\lambda_{\mathrm{b}}}$ where
$\mathbb{S}^{2}\subset\mathbb{R}^{3}$ represents Earth’s spherical surface. In
this model, we assume that each device connects to its nearest BS. Since both
the devices and BS are located on the Earth’s surface, the distances between
the users and their serving BS are very small, in comparison to the Earth’s
circumference, therefore we can approximate the Earth surface as flat and
infinite from the user perspective. As a result, the contact-distance is
simplified to the well-known 2-dimensional PPP model, where the the PDF of the
contact distance $R_{\mathrm{o}}$ is given by [7],

$f_{R_{\mathrm{o}}}(r)=2\pi\lambda_{\mathrm{b}}r\exp(-\pi\lambda_{\mathrm{b}}r^{2}),$
(4)

where $r$ is the distance between the user and its serving BS. Moreover, the
users are assumed to be randomly and homogeneously dispersed on the Earth’s
surface. Hence, we represent their locations with a third independent PPP on
Earth’s surface, i.e. $\Phi_{\mathrm{d}}$ on $\mathbb{S}^{2}$ with density
$\lambda_{\mathrm{d}}$. In this proposed framework we assume that both the
satellite constellation and terrestrial BS listen to the frames transmitted by
users. As such, when a user transmits a frame it is suffice for either of the
networks to receive it for the transmission to be successful. Fig 2
illustrates a user-centric representation of the hybrid network showing the
main geometric parameters.

###  II-B Channel Model

We assume that users are equipped with hemispheric-gain antennas (pointed
upward) such that the uplink EIRP is uniform from the satellites’ perspective,
denoted as $P$. We have two channel models in a hybrid network; (i) ground-to-
satellite and (ii) user-to-terrestrial BS. For the Satellite-to-Ground
channel, we follow the empirical model obtained from our previous work in [6],
in which the received power at the satellite is modeled as,

$S_{\mathrm{s}}=P~{}l~{}\zeta,$ (5)

where $l$ is the free-space path-gain obtained as follows,

$l(\varphi)=\frac{l_{\mathrm{o}}l_{\mathrm{air}}}{d^{2}}=\frac{l_{\mathrm{o}}l_{\mathrm{air}}}{R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}^{2}+(R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}+h)^{2}-2R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}(R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}+h)\cos\varphi},$
(6)

where $l_{\mathrm{air}}$ is the air absorption attenuation caused by the
resonance of gas and water vapor between Earth’s surface and the satellite,
$l_{\mathrm{o}}$ is the path-gain constant ${l_{\mathrm{o}}=c^{2}/(4\pi
f)^{2}}$ where $c$ is the speed of light and $f$ is the center carrier
frequency. $\zeta$ in (5) is the excess path-gain, i.e. the reciprocal of the
excess path-loss. The distribution of the excess path-gain follows a mixed
Gaussian distribution as follows [6],

$\zeta\mathrm{[dB]}\sim
p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}\mathcal{N}(-\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}},\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}^{2})+p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}\mathcal{N}(-\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}},\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}^{2}),$
(7)

where
$\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}},\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}},\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}},\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}$
depend on the propagation environment,
$p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}$ and
$p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}$ are the line-of-sight
(LoS) and non-line-of-sight (NLoS) probabilities respectively. The LoS
probability is derived as a function of the contact angle as follows,

$p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}=\exp(-\beta\cot\theta)=\exp\left(-\frac{\beta\sin\varphi}{\cos\varphi-\alpha}\right),$
(8)

where $\theta$ is the ground user elevation angle and $\beta$ is a parameter
that depends on the propagation environment. On the other hand, for the
terrestrial channel model we utilize the common log-distance path-loss model,
where the received power at the terrestrial BS is given as follows,

$S_{\mathrm{b}}=Pg~{}b~{}l_{\mathrm{o}}r^{-a},$ (9)

where $a$ is the path-loss exponent, $b$ is a model constant, and $g$ is the
fading process assumed to follow a Rayleigh model.

##  III Uplink Interference

Wireless communications rely on the reuse of spectral resources to achieve
better spectral efficiency, thus networks utilize well-designed _access
systems_ that enable the coexistence of different transmissions. Such access
systems endeavor to mitigate the impact of co-channel interference by
intelligently scheduling spectral resources. We capture the effect of access
systems in mitigating the interference via a factor ${\kappa\in[0,1]}$. This
factor acts on the total interference power, where a value of ${\kappa=0}$
indicates an ideal system, i.e. full interference mitigation, and a value of
${\kappa=1}$ indicates the worst-case scenario, i.e. no interference
mitigation. For a generic network, the effect of the access system in the
satellite link is denoted as, $\kappa_{\mathrm{s}}$ and that of the
terrestrial link is denoted as, $\kappa_{\mathrm{b}}$, where two types of
interference are formulated in a hybrid system; (i) at the satellite and (ii)
at the BS receivers. At a given time instance, a satellite can only serve a
certain region underneath it, named the satellite footprint as depicted in Fig
3, this region is either confined because of the antenna beamwidth or the
Earth’s curvature. We denoted this spherical cap as $\mathcal{A}$ having an
Earth-centered apex angle of $2\varphi_{\mathrm{m}}$ and radius
$R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}$, where $\varphi_{\mathrm{m}}$ is
the maximum Earth-centered zenith angle which is given by,

$\varphi_{\mathrm{m}}=\begin{cases}\operatorname{asin}\left(\frac{1}{\alpha}\sin\frac{\psi}{2}\right)-\frac{\psi}{2},&\psi<2\operatorname{asin}\alpha\\\
\operatorname{acos}\alpha,&\psi\geq 2\operatorname{asin}\alpha\\\
\end{cases}~{},$ (10)

where $\psi$ is the satellite ideal beamwidth. Fig 3 also illustrates the
effect of the beamwidth on $\varphi_{\mathrm{m}}$. From the satellite
perspective, the total interference power acting on the signal transmitted by
a user $x_{\mathrm{o}}$ is given by,

$I_{\mathrm{s}}=\sum_{x_{i}\in\Phi_{\mathrm{d}}\cap\mathcal{A}\backslash
x_{o}}\kappa_{\mathrm{s}}Pl_{i}\zeta_{i}~{}.$ (11)

When considering the large swath of Earth covered by the satellite, the
practical number of interfering devices converges to the mean value,
$\lambda_{\mathrm{d}}||\mathcal{A}||$ if all devices are assumed to be active.
Accordingly, the random variations in the interference power are small and
thus could be ignored. As a result, it is sufficient to describe the
interference by its average power to understand the performance of the system.
By invoking Campbell’s theorem of sums [7], we obtain the average interference
power as follows,

$\displaystyle\bar{I_{\mathrm{s}}}=\mathbb{E}\left[I_{\mathrm{s}}\right]$
$\displaystyle=\mathbb{E}_{\zeta}\left[2\pi
R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}^{2}D\lambda_{\mathrm{d}}\int_{0}^{\varphi_{\text{m}}}\kappa_{\mathrm{s}}Pl(\varphi)\zeta(\varphi)\sin\varphi\mathrm{d}\varphi\right]$
$\displaystyle\stackrel{{\scriptstyle(\mathrm{a})}}{{=}}2\pi
R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}^{2}D\lambda_{\mathrm{d}}P\int_{0}^{\varphi_{\text{m}}}\kappa_{\mathrm{s}}l(\varphi)\bar{\zeta}(\varphi)\sin\varphi\mathrm{d}\varphi,$
(12)

where the average density of active users is a homogeneously thinned process,
$D\lambda_{\mathrm{d}}$ resulting from the spatial duty cycle denoted as $D$
and step (a) results from the fact that the excess path-gain $\zeta$ is
independent from the geometric process. The term, $2\pi
R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}^{2}\sin\varphi~{}\mathrm{d}\varphi$
in (III), represents the circular strip area on Earth’s spherical surface for
a given zenith angle. Since $\zeta$ is modeled as Gaussian mixture model
(GMM), its mean value is derived as follows,

$\bar{\zeta}=p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}\exp\left(\frac{\rho^{2}\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}^{2}}{2}-\rho\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}\right)+p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}\exp\left(\frac{\rho^{2}\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}^{2}}{2}-\rho\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}\right),$
(13)

where $\rho=\ln 10/10$. On the other hand, from the terrestrial BS
perspective, the aggregated interference power acting on the signal
transmitted by a user $x_{\mathrm{o}}$ is given by,

$I_{\mathrm{b}}=\sum_{x_{i}\in\Phi_{\mathrm{d}}\backslash
x_{o}}\kappa_{\mathrm{b}}Ph_{i}r^{-a}~{}.$ (14)

If the path-loss exponent is relatively large, then the interference is
significantly influenced by variations in the number of nearby devices. As
such we deal with the terrestrial interference as a random variable. The
expression of the interference’s Laplace Transform (LT) is derived based on
[7] as follows,

$\mathcal{L}_{I_{\mathrm{b}}}(s)=\exp\left[-\pi
D\lambda_{\mathrm{d}}\frac{(\kappa_{\mathrm{b}}Pbl_{\mathrm{o}}s)^{\frac{2}{a}}}{\operatorname{sinc}\frac{2}{a}}\right],$
(15)

where the LT of $I$ is defined as
${\mathcal{L}_{I}(s)=\mathbb{E}[\exp(-sI)]}$.

Figure 3: A satellite-centric representation of the beamwidth effect on its
footprint and its maximum Earth-centered zenith angle.

##  IV Hybrid Coverage Probability

In the proposed hybrid network model, it is suffice to receive the
transmission from a user by either the satellite constellation or the
terrestrial BS. Thus, a transmission will only fail if both of these two
conditions are satisfied: (i) the signal-to-interference-and-noise ratio
(SINR) of the satellite link $\gamma_{\mathrm{s}}$ is below the service
threshold $\gamma_{\mathrm{o}}$, and (ii) the SINR of the terrestrial link is
also below $\gamma_{\mathrm{o}}$. Accordingly, the success probability of the
hybrid coverage for a target SINR is formulated as follows,

$\displaystyle\bar{p}_{\mathrm{c}}$
$\displaystyle=1-\mathbb{E}\left[\mathbb{P}(\gamma_{\mathrm{b}}<\gamma_{\mathrm{o}},\gamma_{\mathrm{s}}<\gamma_{\mathrm{o}})\right]$
$\displaystyle\stackrel{{\scriptstyle(\mathrm{a})}}{{=}}1-\mathbb{E}\left[\mathbb{P}(\gamma_{\mathrm{s}}<\gamma_{\mathrm{o}})\right]~{}\mathbb{E}\left[\gamma_{\mathrm{b}}<\gamma_{\mathrm{o}})\right]$
$\displaystyle=1-\left[1-\bar{p}_{\mathrm{s}}(N_{\mathrm{s}})\right]~{}\left[1-\bar{p}_{\mathrm{b}}(\lambda_{\mathrm{b}})\right]~{},$
(16)

where step (a) stems from the fact that both the channel and geometric
processes of the satellite and terrestrial networks are mutually independent.

Figure 4: Average hybrid probability of coverage for variable
$\lambda_{\mathrm{b}}$ showing the impact of the constellation size
$N_{\mathrm{s}}$, comparing analytic and simulation results. Figure 5: Average
hybrid coverage probability for variable $\lambda_{\mathrm{d}}$ showing the
impact of the constellation size $N_{\mathrm{s}}$, comparing analytic and
simulation results.

The SINR of the received signal at the serving satellite is defined as,
$\gamma_{\mathrm{s}}=S_{\mathrm{s}}/(I_{\mathrm{s}}+W_{\mathrm{s}})$, where
$W_{\mathrm{s}}$ is the average noise power observed at the satellite
receiver. The resulting satellite coverage probability is thus derived as,

$\displaystyle\bar{p}_{\mathrm{s}}(N_{\mathrm{s}})$
$\displaystyle=\mathbb{E}_{\varphi_{\mathrm{o}}}\left[\mathbb{P}\left(\zeta>\frac{\gamma_{\mathrm{o}}[\bar{I}_{\mathrm{s}}+W_{\mathrm{s}}]}{Pl(\varphi_{\mathrm{o}})}\right)\right]$
$\displaystyle=\mathbb{E}_{\varphi_{\mathrm{o}}}\left[1-F_{\zeta}\left(\frac{\gamma_{\mathrm{o}}[\bar{I}_{\mathrm{s}}+W_{\mathrm{s}}]}{Pl(\varphi_{\mathrm{o}})}\right)\right]$
$\displaystyle=\int_{0}^{\varphi_{\text{m}}}\left[1-F_{\zeta}\left(\frac{\gamma_{\mathrm{o}}[\bar{I}_{\mathrm{s}}+W_{\mathrm{s}}]}{Pl(\varphi)}\right)\right]f_{\varphi_{\mathrm{o}}}(\varphi)\mathrm{d}\varphi,$
(17)

where $F_{\zeta}(.)$ is the CDF of the mixed GMM path-gain random variable (in
linear form) is obtained as follows,

$\displaystyle F_{\zeta}(x)=\frac{1}{2}+$
$\displaystyle\frac{p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}}{2}\text{erf}\left[\frac{10\log_{10}x+\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}}{\sqrt{2}\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{LoS}$}}}}}\right]$
$\displaystyle+\frac{p_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}}{2}\text{erf}\left[\frac{10\log_{10}x+\mu_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}}{\sqrt{2}\sigma_{\vbox{\hbox{\scalebox{0.5}{$\mathrm{NLoS}$}}}}}\right].$
(18)

In regards to the terrestrial network, the SINR of the received signal at the
serving BS is defined as,
${\gamma_{\mathrm{b}}=S_{\mathrm{b}}/(I_{\mathrm{b}}+W_{\mathrm{b}})}$, where
$W_{\mathrm{b}}$ is the terrestrial BS’s average noise power. The resulting
terrestrial coverage probability is derived as,

$\displaystyle\bar{p}_{\mathrm{b}}$
$\displaystyle(\lambda_{\mathrm{b}})=\mathbb{E}_{R_{\mathrm{o}}}\left[\mathbb{P}\left(g>\frac{\gamma_{\mathrm{o}}r^{a}(I_{\mathrm{b}}+W_{\mathrm{b}})}{Pbl_{\mathrm{o}}}\right)\right]$
$\displaystyle=\mathbb{E}_{R_{\mathrm{o}}}\left[\mathbb{E}_{I_{\mathrm{b}}}\left[\exp\left(-\frac{\gamma_{\mathrm{o}}r^{a}I_{\mathrm{b}}}{Pbl_{\mathrm{o}}}\right)\right]\exp\left(-\frac{\gamma_{\mathrm{o}}r^{a}W_{\mathrm{b}}}{Pbl_{\mathrm{o}}}\right)\right]$
$\displaystyle=\int_{0}^{\infty}\mathcal{L}_{I_{\mathrm{b}}}\left(\frac{\gamma_{\mathrm{o}}r^{a}}{Pbl_{\mathrm{o}}}\right)\exp\left(-\frac{\gamma_{\mathrm{o}}r^{a}W_{\mathrm{b}}}{Pbl_{\mathrm{o}}}\right)f_{R_{\mathrm{o}}}(r)\mathrm{d}r.$
(19)

We depict in Fig. 4 the coverage probability for different design terrestrial
BS density $\lambda_{\mathrm{b}}$. It is clear from the figure that when we
have low BS density, the satellite network takes over and dominates the hybrid
service and the coverage probability converges to that of the satellite
coverage probability. In order to visualize the impact of user’s density
$\lambda_{\mathrm{d}}$ we depict in Fig. 5 the hybrid coverage with different
satellites constellation sizes. It is clear that as the number of users
increases, the coverage probability deteriorates because of the increasing
interference. With the aid of densifying the satellite constellation, the
coverage probability can be improved again since more serving satellites are
becoming closer to the zenith, i.e. smaller average contact angle leads to
enhanced received powers.

Figure 6: Hybrid operating curves for different QoS levels.

##  V Network Design Parameters

A network design involves determining the network infrastructure parameters to
maintain a target QoS, i.e. coverage probability, $\bar{p}_{\mathrm{c}}$, for
a given user density. While a typical terrestrial network requires evaluating
the BS density only, a hybrid network requires evaluating a combination of the
number of satellites as well as the BS density. Thus, a designer has the
choice of different combinations of $N_{\mathrm{s}}$ and
$\lambda_{\mathrm{b}}$ for a given target performance. In Fig. 6 we depict an
example of such combinations, we call the resulting curve by hybrid network
operating curve, where each curve represents a desired QoS or
$\bar{p}_{\mathrm{c}}$. These curves are crucial to inform strategic
deployment and expansion decisions of such networks, for a given BS density,
the minimum number of satellites required to maintain a certain QoS, i.e. a
point on desired curve, is obtained as follows,

$N_{\mathrm{s}}^{\star}=\operatorname*{arg\,min}_{N_{\mathrm{s}}}\left|\bar{p}_{\mathrm{s}}-\frac{\bar{p}_{\mathrm{c}}-\bar{p}_{\mathrm{b}}(\lambda_{\mathrm{b}})}{1-\bar{p}_{\mathrm{b}}(\lambda_{\mathrm{b}})}\right|~{},$
(20)

where $N_{\mathrm{s}}^{\star}$ is the number of satellites required to
maintain a given $\bar{p}_{\mathrm{c}}$ for a fixed $\lambda_{\mathrm{b}}$. On
the other hand, for a fixed number of satellites, we obtain the BS density
that maintains a certain $\bar{p}_{\mathrm{c}}$ as follows,

$\lambda_{\mathrm{b}}^{\star}=\operatorname*{arg\,min}_{\lambda_{\mathrm{b}}}\left|\bar{p}_{\mathrm{b}}-\frac{\bar{p}_{\mathrm{c}}-\bar{p}_{\mathrm{s}}(N_{\mathrm{s}})}{1-\bar{p}_{\mathrm{s}}(N_{\mathrm{s}})}\right|~{}.$
(21)

The hybrid operating curve can be obtained for a QoS level by varying either
$N_{\mathrm{s}}$ or $\lambda_{\mathrm{b}}$ for a fixed $\bar{p}_{\mathrm{c}}$.
Both (20) and (21) can be evaluated by utilizing numerical optimization
methods such as the Quasi-Newton algorithm [8]. The coverage of PPP satellite
constellations behaves as a lower bound or worst-case scenario when compared
to more regular constellations such as Walker-delta (e.g. Starlink) and
Walker-star (e.g. OneWeb). As such, these operating curves can provide
designers with bounds for an appropriate expansion strategy for different
densities of users. This can be seen in Fig. 7 where the operating curves are
compared between the PPP, Walker-delta, and Walker-star constellations for an
example height of 500 km. The number of planes for both Walker constellations
is assigned to $\sqrt{N_{\mathrm{s}}}$ such that the number of satellites per
plane is equal to the number of orbital planes. On the other hand, the
inclination angles for Walker-delta and Walker-star constellations are chosen
as 86.4∘ and 53∘ respectively. These operating curves provide a rapid insight
into the performance of the network such that for an increase in the number of
devices, one can opt to increase $N_{\mathrm{s}}$ only (vertical curve shift),
increase $\lambda_{\mathrm{b}}$ (horizontal curve shift), or increase both
(region between vertical and horizontal shifts) as illustrated by the dashed
lines in Fig. 7.

Figure 7: Hybrid operating curve at a target $\bar{p}_{\mathrm{c}}=0.8$
showing the impact of increasing $\lambda_{\mathrm{d}}$, comparing PPP,
Walker-delta, and Walker-star constellations.

##  VI Notes on Simulation

A Monte Carlo simulator is developed to validate the proposed framework. It
generates three independent random processes as follows; (i) a process that
represents the randomly generated satellite points on a sphere with an
altitude $R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}+h$, (ii) a process that
represents the terrestrial BS, randomly scattered on Earth’s surface, and
(iii) a process that represents the users, also randomly scattered on Earth’s
surface. Note that the simulator assumes the Earth as an ideal sphere with an
average radius of $R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}$. In order to
achieve a uniform distribution of points on the spherical surface, the
coordinates are generated as, $\vartheta\sim\mathcal{U}(0,2\pi)$ and
$\phi\sim\operatorname{asin}[\mathcal{U}(-1,1)]$ where $\mathcal{U}(a,b)$ is
the uniform random variable between $a$ and $b$. Table I summarizes the
parameters used in the simulation. The center frequency is set at 2 GHz since
we follow the channel model parameters obtained based on real measurements in
[6]. Furthermore, shared spectrum IoT is limited to low frequency bands such
as the L and S bands by spectrum regulators [9]. At these low frequencies, we
can assume that the air absorption attenuation is negligible based on the ITU
model [10].

TABLE I: Notations and Symbols Symbol | Value [Unit] | Definition
---|---|---
$R_{\vbox{\hbox{\scalebox{0.5}{$\oplus$}}}}$ | 6371 [km] | Earth’s average radius
$h$ | 500 [km] | Satellite orbit height
$f$ | 2 [GHz] [6] | Center frequency
$l_{\mathrm{air}}$ | 0 [dB] [10] | Air absorption model
$\beta$ | 2.3 [6] | LoS probability parameter
$\mu_{\mathrm{LoS}}$ | 0 [dB] [6] | LoS excess path-loss mean
$\sigma_{\mathrm{LoS}}$ | 2.8 [dB] [6] | LoS excess path-loss standard deviation
$\mu_{\mathrm{NLoS}}$ | 12 [dB] [6] | NLoS excess path-loss mean
$\sigma_{\mathrm{NLoS}}$ | 9 [dB] [6] | NLoS excess path-loss standard deviation
$a$ | 3.68 [11] | Path-loss exponent
$P$ | 23 [dBm] | Device transmit EIRP
$\psi$ | 2$\pi$ (isotropic) | Satellite beamwidth
$\gamma_{\mathrm{o}}$ | -20 [dB] | Target (threshold) SINR
$W_{\mathrm{b}}$ | -117 [dBm] | BS average noise power
$W_{\mathrm{s}}$ | -130 [dBm] | Satellite average noise power
$\kappa_{\mathrm{b}}$ | -20 [dB] | Terrestrial interference mitigation factor
$\kappa_{\mathrm{s}}$ | -20 [dB] | Satellite interference mitigation factor
$D$ | 1% | Transmission duty cycle

##  VII Conclusion

This paper presented an analytical approach for capturing the performance of a
hybrid uplink network composing of terrestrial BS and an overlaying satellite
constellation. The approach utilizes tools from stochastic geometry to capture
the system as three independent random processes. The formulated framework
provides designers with tools to evaluate the needed network infrastructure.
The paper presents the effect of the increase of user density on the required
network parameters to maintain the target quality of service allowing the
projection of the future performance of the network.

## References

  * [1] M. Jia, X. Gu, Q. Guo, W. Xiang, and N. Zhang, “Broadband hybrid satellite-terrestrial communication systems based on cognitive radio toward 5G,” _IEEE Wireless Communications_ , vol. 23, no. 6, pp. 96–106, 2016\.
  * [2] B. Evans, M. Werner, E. Lutz, M. Bousquet, G. E. Corazza, G. Maral, and R. Rumeau, “Integration of satellite and terrestrial systems in future multimedia communications,” _IEEE Wireless Communications_ , vol. 12, no. 5, pp. 72–80, 2005.
  * [3] A. Al-Hourani, “An analytic approach for modeling the coverage performance of dense satellite networks,” _IEEE Wireless Communications Letters_ , vol. 10, no. 4, pp. 897–901, 2021.
  * [4] N. Okati, T. Riihonen, D. Korpi, I. Angervuori, and R. Wichman, “Downlink coverage and rate analysis of low earth orbit satellite constellations using stochastic geometry,” _IEEE Transactions on Communications_ , vol. 68, no. 8, pp. 5120–5134, 2020.
  * [5] A. Yastrebova, I. Angervuori, N. Okati, M. Vehkaperä, M. Höyhtyä, R. Wichman, and T. Riihonen, “Theoretical and simulation-based analysis of terrestrial interference to LEO satellite uplinks,” in _GLOBECOM 2020 - 2020 IEEE Global Communications Conference_ , 2020, pp. 1–6.
  * [6] A. Al-Hourani and I. Guvenc, “On modeling satellite-to-ground path-loss in urban environments,” _IEEE Communications Letters_ , vol. 25, no. 3, pp. 696–700, 2021.
  * [7] M. Haenggi, _Stochastic Geometry for Wireless Networks_. Cambridge University Press, 2013.
  * [8] D. F. Shanno, “Conditioning of quasi-newton methods for function minimization,” _Mathematics of Computation_ , vol. 24, no. 111, pp. 647–656, 1970.
  * [9] B. Al Homssi, A. Al-Hourani, Z. Krusevac, and W. S. T. Rowe, “Machine learning framework for sensing and modeling interference in IoT frequency bands,” _IEEE Internet of Things Journal_ , vol. 8, no. 6, pp. 4461–4471, 2021.
  * [10] “Propagation Data Required for the Design of Earth-Space Land Mobile Telecommunication Systems,” International Telecommunication Union (ITU), P Series document ITU-R P.681-6, Tech. Rep., 2019.
  * [11] B. Al Homssi, A. Al-Hourani, S. Chandrasekharan, K. M. Gomez, and S. Kandeepan, “On the bound of energy consumption in cellular IoT networks,” _IEEE Transactions on Green Communications and Networking_ , vol. 4, no. 2, pp. 355–364, 2020.


Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space and $X$ be a
random variable in this space. For a sub-sigma algebra
$\mathcal{F}^{\prime}\subseteq\mathcal{F}$ the _conditional expectation of $X$
given $\mathcal{F}^{\prime}$_ is an $\mathcal{F}^{\prime}$-measurable random
variable denoted $\mathbb{E}[X\mid\mathcal{F}^{\prime}]$ such that for every
set $A\in\mathcal{F}^{\prime}$ it holds
$\mathbb{E}[X\cdot\mathbb{I}_{A}]=\mathbb{E}[\mathbb{E}[X\mid\mathcal{F}^{\prime}]\cdot\mathbb{I}_{A}]$,
where $\mathbb{I}_{A}$ is the indicator function of the set $A$. There may
generally be zero or multiple $\mathcal{F}^{\prime}$-measurable variables
satisfying the defining condition of conditional expectation, but if at least
one exists, all of the others differ from it only a set of zero probability.
Hence, when at least one such random variable exists, any of them can be
picked as $\mathbb{E}[X\mid\mathcal{F}^{\prime}]$.

If $B\in\mathcal{F}$ is an event of positive probability, the conditional
expectation of $X$ given $B$ is the expectation of $X$ w.r.t. the probability
measure $\mathbb{P}[\leavevmode\nobreak\ \cdot\mid
B]=\frac{\mathbb{P}[\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ \cap
B]}{\mathbb{P}[B]}$.

###### Definition C.1.

Let $(Y_{i})_{i=0}^{\infty}$ be a (1-dimensional) stochastic process adapted
to some filtration $(\mathcal{F}_{i})_{i=0}^{\infty}$ such that
$\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}]$ exists for every $i\geq 0$. We say
that the process is a _supermartingale_ if for every $i\geq 0$ it holds
$\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}]\leq Y_{i}$. We call the process a
_submartingale_ if $\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}]\geq Y_{i}$ for
every $i\geq 0$.

###### Definition C.2 (Stopping time).

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space and
$(\mathcal{F}_{i})_{i=0}^{\infty}$ a filtration. A _stopping time_ is a random
variable $T\colon\Omega\rightarrow\mathbb{N}\cup\\{\infty\\}$ s.t.
$\\{\omega\in\Omega\mid T(\omega)\leq t\\}\in\mathcal{F}_{t}$ for any
$t\in\mathbb{N}$.

###### Theorem C.3 (Optional stopping theorem, OST).

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space. Next, let
$(Y_{i})_{i=0}^{\infty}$ be a 1-dimensional stochastic process adapted to some
filtration $(\mathcal{F}_{i})_{i=0}^{\infty}$ of $\mathcal{F}$, and $T$ be a
stopping time w.r.t. the same filtration $(\mathcal{F}_{i})_{i=0}^{\infty}$.
Assume that $E[|Y_{i}|]<\infty$ for all $i\geq 0$ and that the above objects
satisfy one of the following conditions:

  * (C1’)

There exists a constant $c$ such that $T\leq c$ with probability 1 (i.e., the
stopping time is almost-surely bounded).

  * (C2’)

There exists a constant $c$ such that for each $t\in\mathbb{N}$ and each
$\omega\in\Omega$ it holds $|Y_{\min\\{t,T(\omega)\\}}(\omega)|\leq c$ (i.e.,
the process is bounded from both below and above up until the point of
stopping).

  * (C3’)

$\mathbb{E}[T]<\infty$, $\mathbb{E}[|Y_{0}|]<\infty$, and there exists a
constant $c$ such that for every $t\in\mathbb{N}$ it holds
$\mathbb{E}[|Y_{t+1}-Y_{t}|\mid\mathcal{F}_{t}]\leq c$ (i.e., the expected
one-step change of the process is uniformly bounded over its evolution, even
if conditioned by the whole past history of the process).

Then, $\mathbb{E}[Y_{T}]$ is well-defined, and moreover
$\mathbb{E}[Y_{T}]\leq\mathbb{E}[Y_{0}]$ if $(Y_{i})_{i=0}^{\infty}$ is a
supermartingale and $\mathbb{E}[Y_{T}]\geq\mathbb{E}[Y_{0}]$ if
$(Y_{i})_{i=0}^{\infty}$ is a submartingale. In other words, the expected
value of, say supermartingale, at the point of stopping is bounded from above
by its mean initial value, and dually for submartingales.

We conclude this section by restating the Extended optional stopping theorem
from (Wang et al., 2019).

###### Theorem C.4 (Extended OST, (Wang et al., 2019)).

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space. Next, let
$(Y_{i})_{i=0}^{\infty}$ be a 1-dimensional stochastic process adapted to some
filtration $(\mathcal{F}_{i})_{i=0}^{\infty}$ of $\mathcal{F}$, and $T$ be a
stopping time w.r.t. the same filtration $(\mathcal{F}_{i})_{i=0}^{\infty}$.
Assume that the above objects satisfy the following condition:

  * (C4’)

There exist real numbers $M,c_{1},c_{2},d$ such that (i) for all sufficiently
large $n\in\mathbb{N}$ it holds $\mathbb{P}(T>n)\leq c_{1}\cdot e^{-c_{2}\cdot
n}$; and (ii) for all $t\in\mathbb{N}$ it holds $|Y_{n+1}-Y_{n}|\leq M\cdot
n^{d}$.

Then, $\mathbb{E}[Y_{T}]$ is well-defined, $\mathbb{E}[|Y_{i}|]<\infty$ for
every $i\geq 0$ and moreover, $\mathbb{E}[Y_{T}]\leq\mathbb{E}[Y_{0}]$ if
$(Y_{i})_{i=0}^{\infty}$ is a supermartingale and
$\mathbb{E}[Y_{T}]\geq\mathbb{E}[Y_{0}]$ if $(Y_{i})_{i=0}^{\infty}$ is a
submartingale.

##  Appendix D Proof of Theorem 5.4

See 5.4

###### Proof.

We present the proof for $L_{f}$, the proof for $U_{f}$ is analogous.

Recall that $Z_{n}$ denotes the $n$-th state along a run of $\mathcal{C}$ and
$\mathbf{X}_{n}$ denotes the $n$-th valuation encountered along a run. Let us
define a stochastic process $Y=(Y_{n})_{n=0}^{\infty}$ by putting
$Y_{n}\mathrel{\mathop{\ordinarycolon}}=L_{f}(Z_{n})+f(\mathbf{X}_{n}^{\mathit{out}})$.
The process $Y$ is clearly adapted to the canonical filtration
$(\mathcal{F}_{n})_{n=0}^{\infty}$.

First, note that $\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}]$ exists and for any
run $\rho$ is defined as follows: let $\tau=(\ell,\mathit{Pr})$ be the unique
transition enabled in $Z_{i}(\rho)$. Then

(7) $\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}](\rho)=\sum_{\ell^{\prime}\in
L}\mathit{Pr}(\ell^{\prime})\cdot\mathbb{E}[L_{f}(\ell^{\prime},\mathbf{N})+f(\mathbf{N}^{\mathit{out}}))],$

where $\mathbf{N}=\mathit{Next}(\tau,\mathbf{X}_{i}(\rho))$. This function is
well-defined, since the expectation on the right-hand side of (7) always
exists by the OST-soundness assumption. In Section E, we prove that the
function defined in this way indeed satisfies the definition of conditional
expectation.

Next, we prove that the process $Y$ is a submartingale. We can continue from
(7) as follows:

$\displaystyle\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}](\rho)$
$\displaystyle=\sum_{\ell^{\prime}\in
L}\mathit{Pr}(\ell^{\prime})\cdot\mathbb{E}[L_{f}(\ell^{\prime},\mathbf{N})+f(\mathbf{N}^{\mathit{out}}))]$
(by (7)) $\displaystyle\geq L_{f}(Z_{i}(\rho))+f(Z_{i}(\rho))$ (by (4))
$\displaystyle=Y_{i}(\rho)$ (by the def. of $Y_{i}$),

as required.

In what follows, we abbreviate $\mathit{TimeTerm}$ by $T$. Since
$(\mathcal{C},L_{f},f)$ is OST-sound, the submartingale $Y$ satisfies the
assumptions of either the optional stopping theorem or its extended variant.
It follows that

(8) $\displaystyle
L_{f}(\ell_{\mathit{init}},\mathbf{x}_{\mathit{init}})+f(\mathbf{x}_{\mathit{init}}^{\mathit{out}})=\mathbb{E}[Y_{0}]\leq\mathbb{E}[Y_{T}]=\mathbb{E}[f(\mathbf{X}_{T}^{\mathit{out}})]=\mathbb{E}_{\mathbf{x}\sim\mu^{\mathcal{C}}}[f(\mathbf{x}^{\mathit{out}})],$

where the first equality follows from the definition of $Y$, the second from
the (extended) optional stopping theorem, the third from the fact that $L_{f}$
is zero upon termination, and last one from the definition of
$\mu^{\mathcal{C}}$. ∎

##  Appendix E Conditional Expectation for U/LESMs

We argue that

(9) $\mathbb{E}[Y_{i+1}\mid\mathcal{F}_{i}](\rho)=\sum_{\ell^{\prime}\in
L}\mathit{Pr}(\ell^{\prime})\cdot\mathbb{E}[L_{f}(\ell^{\prime},\mathbf{N})+f(\mathbf{N}^{\mathit{out}}))],$

where $\tau$ is the unique transition enabled in state $Z_{i}(\rho)$.

For each $i\geq 0$ we write $Y_{i}=g(Z_{i})$ for a Borel-measurable function
$g$. Note that the right-hand side of (9) can be written in measure-
theoretical terms as

$\int g(s)dP_{Z_{i}(\rho)}(s),$

where $P_{x}$ is the probability measure on states of the program defined by
the transition kernel of the process represented by our program in the source
state $x$.

Hence, our aim is to prove that for any $\mathcal{F}_{i}$-measureable set $A$
it holds

(10) $\int_{A}g(Z_{i+1}(\rho))d\mathbb{P}(\rho)=\int_{A}\left[\int
g(s)dP_{Z_{i}(\rho)}(s)\right]d\mathbb{P}(\rho),$

where $\mathbb{P}$ is the probability measure over the runs of the program.

We first prove the equality for the case when $A$ is an
$\mathcal{F}_{i}$-_cylinder_ , i.e. a set of the form $A=\\{\rho\mid
Z_{1}(\rho)\in S_{1},\ldots,Z_{i}(\rho)\in S_{i}\\}$ for some Borel-measurable
sets of program states $S_{1},\ldots,S_{i}$. In such a case, the right-hand
side in (10) can be rewritten as

$\int_{S_{1}}\cdots\int_{S_{i}}g(s_{i+1})dP_{s_{i}}{(s_{i+1})}\cdots
dP_{s_{0}}(s_{1}),$

which equals the left-hand side of (10) directly by the cylinder construction
of the probability measure $\mathbb{P}$.

Now to prove that (10) holds for any $\mathcal{F}_{i}$-measurable set $A$,
assume first that $g$ is non-negative. By our OST assumption, both integrals
in (10) are finite. Moreover, the integrals are sigma-additive and an integral
over an empty set is zero. Hence, both integrals define a finite measure over
$\mathcal{F}_{i}$ (the measure of set $A$ being the value of the respective
integral when integrating over $A$). As shown in the previous paragraph, these
two measures agree on the generators of $\mathcal{F}_{i}$ (the
$\mathcal{F}_{i}$-cylinders), and the set of these generators is a
$\pi$-system (is closed under finite intersections). Hence, the two measures
are equal on whole $\mathcal{F}_{i}$ (Williams, 1991, Lemma 1.6). Hence, the
integrals are the same for all $A\in\mathcal{F}_{i}$.

For general $g$, we use the standard trick of splitting $g$ into the non-
negative and negative part: $g=g^{+}-g^{-}$, where both $g^{+}$ and $g^{-}$
are non-negative, and hence integrate, on both sides of (10), to the same
value as shown in the previous paragraph (and this value is finite by the
integrability entailed by OST-soundness), irrespective of the choice of $A$.
The equality of both integrals for $g$ then follows from the linearity of
integrals.

##  Appendix F Algorithm for the Similarity Refunation Problem

We now show how our algorithm can be extended for the Similarity refutation
problem.

#### Additional algorithm parameters.

Recall from Section 4 that the Similarity refutation problem is also defined
with respect to a metric $d$ over the output space
$\mathbb{R}^{|V_{\mathit{out}}|}$ and a lower bound $\epsilon>0$ on the
Kantorovich distance that we wish to prove. Our algorithm inputs $d$ and
$\epsilon$ as parameters. We allow any of the following standard metrics:

  * •

$L^{p}$-metric. We allow the $L^{p}$-metric for any natural number
$p\in\mathbb{N}$. This encapsulates the standard $L^{1}$-metric (i.e.
Manhattan metric) and $L^{2}$-metric (i.e. Euclidean metric). Given
$p\in\mathbb{N}$, the $L^{p}$-metric
$d_{p}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{|V_{\mathit{out}}|}\times\mathbb{R}^{|V_{\mathit{out}}|}\rightarrow\mathbb{R}$
is defined via
$d_{p}(x,y)=(\sum_{i=1}^{|V_{\mathit{out}}|}(|x[i]-y[i]|)^{p})^{1/p}$.

  * •

Discrete metric. We also allow $d$ to be the discrete metric, giving rise to
Total Variation distance between output distributions (see Section 3.3).
Recall, the discrete metric
$d_{0}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{|V_{\mathit{out}}|}\times\mathbb{R}^{|V_{\mathit{out}}|}\rightarrow\mathbb{R}$
is defined via $d_{0}(x,y)=0$ if $x=y$ and $d_{0}(x,y)=1$ if $x\neq y$.

  * •

Uniform metric. Finally, we allow the $L^{\infty}$-metric (i.e. uniform
metric). The uniform metric
$d_{\infty}\mathrel{\mathop{\ordinarycolon}}\mathbb{R}^{|V_{\mathit{out}}|}\times\mathbb{R}^{|V_{\mathit{out}}|}\rightarrow\mathbb{R}$
is defined via $d_{\infty}(x,y)=\max_{1\leq i\leq n}|x[y]-y[i]|$.

#### Algorithm.

Recall from Theorem 5.5 that $f$, $U^{1}_{f}$ and $L^{2}_{f}$ also yield a
lower bound on Kantorovich distance between two output distributions, if we in
addition constrain $f$ to be $1$-Lipschitz continuous over reachable output
sets $I^{1}(\ell_{\mathit{out}}^{1})$ and $I^{2}(\ell_{\mathit{out}}^{2})$ of
the two pCFGs. Hence, our algorithm for the Similarity refutation problem
proceeds analogously as in Section 6.1, with the only difference being that it
collects two additional constraints in Step 2 of the algorithm:

  1. (1)

Lower bound on Kantorovich distance. The algorithm collects the similarity
refutation constraint, according to Theorem 5.5:

$L_{f}(\ell_{\mathit{init}}^{2},\mathbf{x}_{\mathit{init}}^{2})+f((\mathbf{x}_{\mathit{init}}^{2})^{\mathit{out}}-U_{f}(\ell_{\mathit{init}}^{1},\mathbf{x}_{\mathit{init}}^{1})-f((\mathbf{x}_{\mathit{init}}^{1})^{\mathit{out}})\geq\epsilon$

  2. (2)

$1$-Lipschitz continuity. Depending on the choice of the metric $d$ over the
output space, the algorithm impose the $1$-Lipschitz continuity constraint as
follows:

     * •

$1$-Lipschitz continuity w.r.t. $d_{p}$. Suppose that metric $d$ is the
$L^{p}$-metric $d_{p}$ for some $p\in\mathbb{N}$. We enforce the following
constraint for each pCFG $\mathcal{C}_{i}$, $i\in\\{1,2\\}$:

$\begin{split}\forall\mathbf{x},\mathbf{y},\mathbf{a}\in\mathbb{R}^{|V_{\mathit{out}}|}.\,&\mathbf{x},\mathbf{y}\models
I^{i}(\ell_{\mathit{out}}^{i})\land\bigwedge_{j=1}^{|V_{\mathit{out}}|}\Big{(}\mathbf{x}[j]-\mathbf{y}[j]\leq\mathbf{a}[j]\land\mathbf{y}[j]-\mathbf{x}[j]\leq\mathbf{a}[j]\Big{)}\\\
&\Longrightarrow(f(\mathbf{x})-f(\mathbf{y}))^{p}\leq\sum_{j=1}^{|V_{\mathit{out}}|}\mathbf{a}[j]^{p}.\end{split}$

The inequality on the right-hand-side is imposed for all
$\mathbf{a}[j]\geq|\mathbf{x}[j]-\mathbf{y}[j]|$, which is equivalent to
simply imposing it for $\mathbf{a}[j]=|\mathbf{x}[j]-\mathbf{y}[j]|$, giving
rise to a sound and complete encoding of the $1$-Lipschitz continuity.

     * •

$1$-Lipschitz continuity w.r.t. $d_{0}$. The algorithm collects the following
constraints on $f$ to be $1$-Lipschitz continuous over outputs of the pCFGs.
We enforce the following constraint for each pCFG $\mathcal{C}_{i}$,
$i\in\\{1,2\\}$:

$\forall\mathbf{x},\mathbf{y}\in\mathbb{R}^{|V_{\mathit{out}}|}.\,\mathbf{x},\mathbf{y}\models
I^{i}(\ell_{\mathit{out}}^{i})\Rightarrow f(\mathbf{x})-f(\mathbf{y})\leq 1.$

Indeed, a function $f$ is $1$-Lipschitz continuous if for any distinct
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{|V_{\mathit{out}}|}$ the difference in
the values of $f$ is at most $1$, since $d_{0}(\mathbf{x},\mathbf{y})=1$ for
$\mathbf{x}\neq\mathbf{y}$. This gives rise to a sound and complete encoding
of the $1$-Lipschitz continuity with respect to the discrete metric.

     * •

$1$-Lipschitz continuity w.r.t. $d_{\infty}$. The algorithm collects the
following constraints on $f$ to be $1$-Lipschitz continuous over outputs of
the pCFGs. We enforce the following constraint for each pCFG
$\mathcal{C}_{i}$, $i\in\\{1,2\\}$:

$\begin{split}\forall\mathbf{x},\mathbf{y},\mathbf{a}\in\mathbb{R}^{|V_{\mathit{out}}|}.\,\forall
A\in\mathbb{R}.\,&\mathbf{x},\mathbf{y}\models
I^{i}(\ell_{\mathit{out}}^{i})\land\bigwedge_{j=1}^{|V_{\mathit{out}}|}\Big{(}\mathbf{x}[j]-\mathbf{y}[j]\leq\mathbf{a}[j]\land\mathbf{y}[j]-\mathbf{x}[j]\leq\mathbf{a}[j]\Big{)}\\\
&\bigwedge_{j=1}^{|V_{\mathit{out}}|}\Big{(}\mathbf{a}[j]\leq
A\Big{)}\Longrightarrow f(\mathbf{x})-f(\mathbf{y})\leq A.\end{split}$

The above constraint is a sound and complete encoding of the fact that $f$ is
$1$-Lipschitz continuous over $I^{i}(\ell_{\mathit{out}}^{i})$ for each
$i\in\\{1,2\\}$. On the left-hand-side of the entailment, we use the component
$\mathbf{a}[j]$ for each $1\leq j\leq|V_{\mathit{out}}|$ to bound from above
the absolute difference $|\mathbf{x}[j]-\mathbf{y}[j]|$. Moreover, we use $A$
to bound from above the maximum absolute difference. Thus, the inequality on
the right-hand-side is imposed for all
$A\geq\mathbf{a}[j]\geq|\mathbf{x}[j]-\mathbf{y}[j]|$, which is equivalent to
simply imposing it for $A=\mathbf{a}[j]=|\mathbf{x}[j]-\mathbf{y}[j]|$, giving
rise to a sound and complete encoding of the $1$-Lipschitz continuity with
respect to the uniform metric.

As in Section 6.1, the lower-bound constraint is a linear and purely
existentially quantified constraint over the symbolic template variables since
$\mathbf{x}_{\mathit{init}}^{1}$ and $\mathbf{x}_{\mathit{init}}^{2}$ are
fixed. On the other hand, the $1$-Lipschitz continuity constraints are of the
same form as in eq. (6). Hence, we may proceed analogously as in Steps 3 and 4
in Section 6.1 to translate the collected constraints into an LP instance and
reduce the synthesis to LP solving. The algorithm returns “Not
$\epsilon$-output close” and outputs the computed $f$, $U^{1}_{f}$ and
$L^{2}_{f}$ if the LP is successfully solved, or returns “Unknown” otherwise.

The proof of the following theorem can be found in Section G in the
supplementary material.

###### Theorem F.1 (Correctness of Similarity Refutation).

Suppose that the algorithm outputs “Not $\epsilon$-output close”. Then
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are indeed not $\delta$\- output
close, and $U^{1}_{f}$ and $L^{2}_{f}$ are valid UESM and LESM for $f$,
respectively.

#### Optimization of the Kantorovich distance.

Finally, we note that our algorithm for the Similarity refutation problem
reduces the synthesis of $f$, $U^{1}_{f}$ and $L^{2}_{f}$ to an LP instance
without the optimization objective. Thus, our method can also optimize the
lower bound on the Kantorovich distance by treating $\epsilon$ as a variable
and adding the optimization objective to maximize $\epsilon$.

##  Appendix G Soundness Proofs for Algorithms

See 6.1

###### Proof.

Suppose that the algorithm outputs “Not output-equivalent” and that it
computes a function $f$ over outputs, a state function $U^{1}_{f}$ in
$\mathcal{C}_{1}$ and a state function $L^{2}_{f}$ in $\mathcal{C}_{2}$. We
show that $U^{1}_{f}$ is an UESM for $f$ in $\mathcal{C}_{1}$ and that
$L^{2}_{f}$ is an LESM for $f$ in $\mathcal{C}_{2}$ which together prove that
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are not output-equivalent.

Since the algorithm outputs “Not output-equivalent”, we must have that $f$,
$U^{1}_{f}$ and $L^{2}_{f}$ computed by the algorithm provide a part of the
solution to the system of constraints in Step 3. Thus, it follows by the
correctness of reduction to an LP instance that was established in (Chatterjee
et al., 2016; Asadi et al., 2021; Wang et al., 2019; Zikelic et al., 2022)
that they also provide a solution to the system of constraints collected by
the algorithm in Step 2. But the constraints collected in Step 2 impose the
defining conditions of UESMs as in Definition 5.1, the defining conditions of
LESMs as in Definition 5.2 and OST-soundness conditions as in Definition 5.3.
Note that the absolute value of the sum of the U/LESM and $f$ is finite as
required by OST soundness, since the U/LESMs and $f$ are defined via
polynomial expressions and all program variables have all moments finite at
each step of the program execution. (The latter property follows by a
straightforward induction using the fact that the programs use polynomial
updates and only sample from distributions that have all moments finite.)

Hence, any solution to the system of constraints collected in Step 2 of the
algorithm gives rise to $U^{1}_{f}$ and $L^{2}_{f}$ that are valid UESM and
LESM for $f$, as wanted. Furthermore, by the equivalence refutation constraint
that is also collected in Step 2 and by Theorem 5.5, it follows that whenever
a solution to the system of constraints in Step 2 exists, the two pCFGs are
not output-equivalent. This concludes the proof. ∎

See F.1

###### Proof.

Suppose that the algorithm outputs “Not $\epsilon$-output close” and that it
computes a function $f$ over outputs, a state function $U^{1}_{f}$ in
$\mathcal{C}_{1}$ and a state function $L^{2}_{f}$ in $\mathcal{C}_{2}$. We
show that $U^{1}_{f}$ is an UESM for $f$ in $\mathcal{C}_{1}$ and that
$L^{2}_{f}$ is an LESM for $f$ in $\mathcal{C}_{2}$ which together prove that
$\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are not $\epsilon$-output close.

Since the algorithm outputs “Not $\epsilon$-output close”, we must have that
$f$, $U^{1}_{f}$ and $L^{2}_{f}$ computed by the algorithm provide a part of
the solution to the system of constraints in Step 3. Thus, it follows by the
correctness of reduction to an LP instance that was established in (Chatterjee
et al., 2016; Asadi et al., 2021; Wang et al., 2019; Zikelic et al., 2022)
that they also provide a solution to the system of constraints collected by
the algorithm in Step 2. But the constraints collected in Step 2 impose the
defining conditions of UESMs as in Definition 5.1, the defining conditions of
LESMs as in Definition 5.2 and OST-soundness conditions as in Definition 5.3.
Hence, any solution to the system of constraints collected in Step 2 of the
algorithm gives rise to $U^{1}_{f}$ and $L^{2}_{f}$ that are valid UESM and
LESM for $f$, as wanted. Furthermore, by the similarity refutation constraint
that is also collected in Step 2 and by Theorem 5.5, it follows that whenever
a solution to the system of constraints in Step 2 exists, the two pCFGs are
not $\epsilon$-output close. This concludes the proof. ∎

##  Appendix H Experimental Comparison of OST Conditions

| Name | C1/C4 | C2 | C3
---|---|---|---|---
| Eq. Ref. | T.(s) | Dis. | T.(s) | Eq. Ref | T.(s) | Dis. | T.(s) | Eq. Ref | T.(s) | Dis. | T.(s)
4cm Benchmarks From
(Wang et al., 2019) | Simple Example | ✓ | 0.30 | 266.667 | 0.25 | ✓ | 1.75 | 0.583 | 32.18 | ✓ | 0.46 | 266.667 | 0.48
Nested Loop | ✓ | 0.31 | 50.0 | 0.33 | TO | - | TO | - | TO | - | TO | -
Random Walk | ✓ | 0.23 | 9.0 | 0.22 | TO | - | TO | - | ✓ | 1.13 | 11.25 | 0.59
Goods Discount | ✓ | 0.35 | 0.008 | 0.56 | ✓ | 0.52 | 0.008 | 0.91 | ✓ | 1.51 | 0.008 | 0.82
Pollutant Disposal | ✓ | 0.44 | 0.026 | 0.46 | TO | - | TO | - | TO | - | TO | -
2D Robot | ✓ | 13.90 | - | - | TO | - | TO | - | TO | - | TO | -
Bitcoin Mining | ✓ | 0.25 | 0.05 | 0.22 | ✓ | 0.35 | 0.05 | 0.24 | ✓ | 0.36 | 0.05 | 0.3
Bitcoin Mining Pool | ✓ | 129.98 | 122761.25 | 131.06 | TO | - | TO | - | ✓ | 163.00 | TO | -
Species Fight | ✓ | 0.90 | - | - | TO | - | TO | - | TO | - | TO | -
4cm Benchmarks From
(Kura et al., 2019) | coupon_collector | ✓ | 1.10 | 0.5 | 1.41 | TO | - | TO | - | ✓ | 1.69 | 0.5 | 1.80
coupon_collector4 | ✓ | 70.96 | - | - | TO | - | TO | - | ✓ | 210.63 | TO | -
random_walk_1d_intvalued | ✓ | 0.32 | 1.2 | 0.38 | TO | - | TO | - | ✓ | 0.98 | 2.4 | 1.85
random_walk_1d realvalued | ✓ | 0.27 | 3.841 | 0.50 | TO | - | TO | - | TO | - | TO | -
random_walk_1d_adversary | ✓ | 0.38 | 0.768 | 0.40 | TO | - | TO | - | TO | - | TO | -
random_walk_2d_demonic | ✓ | 0.58 | 0.668 | 0.58 | TO | - | TO | - | TO | - | TO | -
random_walk_2d_variant | ✓ | 0.75 | 0.501 | 0.76 | TO | - | TO | - | TO | - | TO | -

Table 2. Comparison of Different OST conditions applied to the first benchmark
set

Table H shows a comparison between performance of different OST conditions
when applied for refuting equivalence/similarity of our benchmarks. As
expected, (C2) and (C3) are more restrictive, therefore fewer benchmarks could
be refuted by them. Moreover, even when with (C2) or (C3) our tool
successfully refutes equivalence, they take more time than (C1)/(C4). This
shows that although (C2) and (C3) can be applied to a wider range of programs,
whenever (C1) or (C4) are applicable, it is more efficient to use the latter.

Note that all of our benchmarks satisfy the assumptions of (C4) while some
also satisfy (C1). Moreover, both (C1) and (C4) do not impose any constraints
on the generated ESMs, therefore they are presented in the same column in
table H.