# Uncertainty Aware ML-based surrogate models for particle accelerators: A
Study at the Fermilab Booster Accelerator Complex

Malachi Schram *  0<EMAIL_ADDRESS>Kishansingh Rajput *  0 Thomas Jefferson
National Accelerator Laboratory, Newport News, VA 23606, USA  Karthik Somayaji
NS *  0 Peng Li *  0 University of California Department of Electrical and
Computer Engineering, Santa Barbara, CA, 93106, USA  Jason St. John *  0 Fermi
National Accelerator Laboratory, Batavia, Illinois 60510, USA  Himanshu Sharma
*  0 Pacific Northwest National Laboratory, Richland, WA 99354, USA

###### Abstract

Standard deep learning methods, such as Ensemble Models, Bayesian Neural
Networks and Quantile Regression Models provide estimates to prediction
uncertainties for data-driven deep learning models. However, they can be
limited in their applications due to their heavy memory, inference cost, and
ability to properly capture out-of-distribution uncertainties. Additionally,
some of these models require post-training calibration which limits their
ability to be used for continuous learning applications. In this paper, we
present a new approach to provide prediction with calibrated uncertainties
that includes out-of-distribution contributions and compare it to standard
methods on the Fermi National Accelerator Laboratory (FNAL) Booster
accelerator complex.

##  I Introduction

Particle accelerators are complex multi-system machines that include a large
number of variables with non-linear dynamics. To date, accelerator control
systems are manually optimized by experts that are guided by physical
principles whenever possible. Developing high-dimensional, physics-based
models that account for multiple time scales is extremely challenging.
Although, there are accelerator beam models with impressive and improving
precision [1], building fully comprehensive Monte Carlo-based models of an
entire facility is challenging, if not intractable. Data-driven methods, such
as deep neural networks (DNNs), are well suited to capture the dynamics of
these non-linear complex systems. These surrogate models can then be used to
develop new AI-based control systems provided they can inform the optimization
algorithm on how reliable the predictions are.

The recent development of DNNs [2, 3, 4] has proven itself useful for complex
control problems [5, 6, 7, 8]. The use of machine learning for particle
accelerator applications has grown in recent years to include, but is not
limited to, diagnostics [9, 10, 11, 12, 13, 14], anomaly
detection/forecasting/classification [15, 16, 17, 18, 19], and controls [20,
21, 22, 23]. Although these studies have shown some impressive results, the
majority do not include any uncertainty quantification (UQ) to compliment
their predictions. Unfortunately, the use of DNNs for online safety-critical
applications remains limited due to issues such as model explainability, in-
domain and out-of-domain prediction and uncertainties, and uncertainty
calibrations.

In recent years, there has been an increasing amount of effort on estimating
uncertainties in DNNs. A prediction’s uncertainty can be separated into the
model’s intrinsic uncertainty (model uncertainty) and the uncertainty caused
by the data (data uncertainty). The model uncertainty is typically reducible,
within limits, by improving the model architecture and hyper-parameters,
however, the data uncertainty is irreducible. Additionally, uncertainty
estimation originating from Out-Of-Distribution (ODD) samples is critical for
a number of applications, such as using DNNs as a proxy to model a dynamical
system used for system control and/or optimization. A deep learning method
that provides predictive uncertainty is not sufficient for safe decision-
making; a deep learning method with properly calibrated uncertainty is
required.

Recent studies that include data-driven, Machine Learning (ML) based surrogate
models have started to include UQ in their models, such as modeling the FNAL
Booster accelerator complex for a reinforcement learning application [21] and
on uncertainty aware anomalies prediction [15]. Additionally, a recently
published study compared the use of Bayesian Neural Networks (BNN) and
ensemble methods for particle accelerator applications [24]. In this paper, we
compare three different methods to estimate data-related uncertainties for DNN
models as it applies to modeling the FNAL Booster Accelerator Complex. In
Section II, we briefly describe the Fermilab Booster accelerator complex and
the data used for training the DNN models in the context of a control
optimization problem. In Section III, we introduce three methods that estimate
uncertainty quantification for DNNs. In Section IV, we present the performance
of each method for in-distribution and out-of-distribution scenarios. Finally,
we conclude with a summary of our results in Section V

##  II Fermilab Booster Accelerator and Complex

At 15 Hz the Fermilab Booster rapid-cycling synchrotron accelerates each
injected batch of 400 MeV protons to 8 GeV and resets to receive the next
injection. See [20] for detailed discussion. A central component of the
Booster cycle is the Gradient Magnet Power Supply (GMPS), which provides the
synchronously rising and falling electrical current to this circular
synchrotron’s main bending magnets, tracking the energy (and therefore the
proton beam’s magnetic stiffness) upward to extraction before returning to the
injection state.

The throughput efficiency of the synchrotron is sensitive to unwanted
perturbations of the GMPS current, causing the beam’s trajectory to deviate
from the desired path and scrape on apertures. Such perturbations are
understood to be induced by the power supplies of other nearby synchrotrons on
their own cycles in the accelerator complex, temperature variations, 60 Hz
power line frequency meanders, and other accelerator complex nuances. A
proportional–integral–derivative (PID [25, 26]) regulator circuit attempts to
compensate for these perturbations with cycle-by-cyle adjustments to the
minimum of the sinusoidal control signal. Figure 1 shows a schematic overview
of the GMPS control environment. See Figure 2 for a sample distribution of
measured errors for the minimum value of the sinusoidally varying magnet
current.

Machine Learning techniques promise a new avenue for developing more
sophisticated control agents with better overall regulation performance,
allowing a predictive, anticipatory approach not encompassed by the reactive
PID regulation paradigm. As the authors in [20] point out, any ML-based GMPS
regulator which replaces this PID circuit is required to deliver stable, fast
inference times; the data intake, forward inference, and generation of the
resulting control signal must always be complete in less than 66 milliseconds.

Reinforcement Learning (RL) was selected as the approach to train a ML-based
agent to act as the GMPS regulator. RL learns an action policy by training a
model using data describing a system’s states, actions, and the resulting
outcomes. This technique is a natural choice because, once a competent agent
is in operation, its real-world performance can be used to provide updated
model parameters, tracking the slowly changing dynamics of the accelerator
complex. At the outset of RL training, the control agent would be expected to
make egregious mistakes and to learn from them. Thus a surrogate model is
needed, one which captures the dynamics of the Booster GMPS regulator’s
control environment, where the agent can learn from mistakes without risk to
personnel, equipment, or the science program they support around the clock.

Figure 1: Schematic view of the GMPS control environment. The human operator
specifies a target program via the Accelerator Control Network that is
transmitted to the GMPS control board. The FPGA-based control logic utilizes
these settings together with readings from a reference magnet to prescribe a
driving signal to the GMPS. The effect of this prescribed signal on the
bending magnets is measured by an in-series reference magnet, with sampled
readings transmitted back to the GMPS control board. Reference measurements
and prescribed signals may be logged and transmitted over network for later
analysis. Figure 2: Distribution of fractional measured error in the GMPS
current at the minimum value of the magnet current (prefix B: indicates
Booster, I current, MIN miniumum, and ER error), with the non-ML PID regulator
discussed in the text. From [20].

Surrogate model training and testing data were taken with the PID regulation
circuit in operation. Additionally a small amount of data were taken with the
regulation circuit off, or with changes to its coefficients. This choice
maximized the available volume of training and testing data, sampling the
changing response dynamics of the GMPS regulation while minimizing impact on
accelerator operations. The PID regulator circuit’s residual error is
typically only 0.1% of the injection current minimum. Without regulation, the
fitted minimum of the magnetic field may vary from the set point by as much as
a few percent.

Five time series were used to produce the surrogate model, and we use their
names as they are logged by the accelerator control system. (In the
accelerator control system’s data-logging nomenclature, device parameters with
the B: prefix are related to the Booster, whereas device parameters beginning
with I: are related to the Main Injector. “MDAT” denotes the accelerator
(“machine”) data communication broadcast.) B:VIMIN is the compensating
recommendation for the minimum value of the offset-sinusoidal GMPS current,
issued by the GMPS regulator in order to reduce the magnitude of B:IMINER,
which is a measure of the residual error. B:LINFRQ is the measured offset from
the expected 60 Hz line frequency powering the GMPS, in mHz. I:IB and I:MDAT40
provide measurements of the Main Injector bending dipole current at different
points in the circuit and through different communication channels. Among all
candidate device time series analyzed, Main Injector’s power supplies were
shown in [20] to have the highest Granger causality with respect to B:IMINER,
the minimization control objective of this application, whose typical value
distribution is shown in Figure 2.

Table 1 briefly summarizes the parameters of interest used in surrogate
modeling.

Table 1: Description of dataset parameters chosen by experts and later
validated with a causality study. Here,“MI” means Main Injector, “MDAT” means
accelerator (machine) data communication, and device parameters that begin
with B are related to the Booster, whereas device parameters that begin with I
are related to the Main Injector. (Reused with permission from the authors of
[20])

Parameter | Details [Units]
---|---
B:IMINER | Setting-error discrepancy at injection [A]
B:LINFRQ | 60 Hz line frequency deviation [mHz]
B:VIMIN | Compensated minimum GMPS current [A]
I:IB | MI lower bend current [A]
I:MDAT40 | MDAT measured MI current [A]

Data were collected nominally at 15 Hz, and due to clock drift among the
front-ends taking those samples, the data as logged were then time-aligned to
a periodic reference signal. Period 0 (June 3, 2019 to July 11, 2019) was
ended by the annual Summer Shutdown and Maintenance. Period 1 (December 3,
2019 to April 13, 2020) ended when the accelerator operations were suspended
in response to the COVID-19 pandemic. More detail on the collection and
preparation of the data can be found at the Data Descriptor article [27].

##  III Machine Learning Methods

There has been a lot of research in uncertainty quantification for deep
learning models that includes, but not limited to, BNN [28], Deep Quantile
Regression (DQR) [29], and Deep Gaussian Process Approximation models (DGPA)
[30, 31, 32, 33]. In this paper, we do not consider ensemble methods because
these methods require training of multiple models and multiple inferences to
provide an uncertainty estimation, making it computationally expensive, slow,
and memory intensive. For this paper, we implemented models for BNN, DQR, and
DGPA to better understand their performance for in-distribution and out-of-
distribution uncertainty estimation and report the results in the next
section. Our specific effort to develop a new DGPA model is most closely
related to the recent paper on Spectral-normalized Neural Gaussian Process
(SNGP) models [34] for classification. In the following subsections we discuss
the working mechanism of these methods.

###  III.1 Bayesian Neural Network (BNN) Model

Monte-Carlo (MC) Dropout [28] is a commonly used approach to estimate the
prediction uncertainty for deep neural network models. The MC-dropout approach
has shown [28] to overcome the computational challenges of estimating
uncertainty using Bayesian models. The review by Abdar et.al [35] provide a
comprehensive details on estimating uncertainty in deep learning. The initial
application of MC dropout was to over-come over fitting associated with deep
neural networks (DNN) while training.

The MC dropout approach relies on introducing a tunable uncertainty into a
network training process by adding a dropout layer that randomly removes nodes
to the following layer with a set probability ($p$) at each forward pass in
training process. Once the model is trained the probability assigned to the
dropout layer can be calibrated to provide an estimation of the data
uncertainty when using for inference. While training the DNN with dropout, the
units in a layer are randomly dropped to avoid over-fitting. The loss function
of DNN with dropout is described in Equation 1 where, $\widehat{\mathbf{y}}$
is the output of the NN model, $E(\cdot,\cdot)$ is the loss function, $W_{i}$
are the weight matrices of each layer $i=1\ldots L$,
$\mathbf{y}_{i}\in\mathbf{Y}$ are the observed output data corresponding to
input $\mathbf{x}_{i}\in\mathbf{X}$ and $\lambda$ is the weight decay
parameter for the $L_{2}$ regularizing-term in the loss function .

$\mathcal{L}_{\text{dropout
}}:=\frac{1}{N}\sum_{i=1}^{N}E\left(\mathbf{y}_{i},\widehat{\mathbf{y}}_{i}\right)+\lambda\underbrace{\sum_{i=1}^{L}\left(\left\|\mathbf{W}_{i}\right\|_{2}^{2}+\left\|\mathbf{b}_{i}\right\|_{2}^{2}\right)}_{regularizing-
term}$ (1)

The predictive uncertainty with the model is expressed by Equation 2, where
$p\left(y^{*}\mid x^{*},w\right)$ is the model likelihood and $p(w\mid X,Y)$
is the posterior on the model weights. Calculating the posterior on weights
are analytically intractable and hence in the MC-dropout approach this is
modeled by $q(w)$ distribution on the model weights $W_{i}$. The columns of
the weight matrices are randomly set to zero using a Bernoulli distribution
$W_{i}=M_{i}\operatorname{diag}\left(\left[z_{i,j}\right]_{j=1}^{K_{i}}\right)$.
$M_{i}$ is the matrix of variational parameters (readers can find additional
details in [28]). Hence, MC-dropout forward yields a different weight matrix,
with the predictive mean of the model output expressed as Equation 3, where
$T$ is the number of MC iterations. The $\hat{y}^{*}$ is the neural network
output for input $x^{*}$ and $p_{MC}$ is the predictive mean.

$p\left(y^{*}\mid x^{*},X,Y\right)=\int p\left(y^{*}\mid x^{*},w\right)p(w\mid
X,Y)dw$ (2)

$E_{q\left(y^{*}\mid
x^{*}\right)\left(y^{*}\right)}\approx\frac{1}{T}\sum_{t=1}^{T}\hat{y}^{*}\left(x^{*},W_{i}^{t}\right)=p_{MC}\left(y^{*}\mid
x^{*}\right)$ (3)

###  III.2 Deep Quantile Regression (DQR) Model

DQR is a method used to estimate the conditional quantiles of a response
variable distribution that is more robust against outliers in the response
measurements [29]. We define the conditional quantile function for a non-
linear relationship in Equation 4.

$Q_{y}(\tau|x_{t})=G_{\tau}(x_{t},w)$ (4)

Here $G_{\tau}(x_{t},w)$ is a non-linear function that is approximated by a
DNN, $x$ is the input feature vector at time $t$, and $\tau^{th}$ is the
conditional quantile. We develop the DNN model to simultaneously learn
predictions based on a set of defined quantiles. For each defined quantile,
the prediction $y_{i}^{p}$ and outcome $y_{i}$, the regression loss for a
$\tau^{th}$ quantile is given in Equation 5:

$\mathcal{L}(\mathbf{y_{i}^{p}},\mathbf{y_{i}})=\max[\tau(\mathbf{y}_{i}-\mathbf{y}_{i}^{p}),(\tau-1)(\mathbf{y}_{i}-\mathbf{y}_{i}^{p})]$
(5)

As such, DQR provides a comprehensive statistical model that captures non-
linear relationships by providing conditional quantiles along with the median,
in contrast to traditional regression methods. In this paper, we implemented a
deep learning model with convolutional layers, similar to the model
architecture used for the other methods discussed in this paper, however, the
output layers mapping to a dedicated quantile value.

###  III.3 Deep Gaussian Process Approximate (DGPA) Model

Gaussian Process (GP) models uses a kernel function to transform the input
data into some higher dimensional representation. The function utilizes the
point-to-point distance between samples in the new representation to produce
predictions. With this property it’s intrinsically distance aware and can
detect OOD samples based on the distance from training distribution.
Unfortunately, GP models do not scale with large data sets and large feature
dimensions. As such, using GP on high dimensional data usually requires either
dimension reduction, feature extraction or some other form of approximation.
In contrast, DNNs can be very expressive and readily applied to problems with
large data sets and high dimensional feature space. Unfortunately,
deterministic DNNs can make predictions on samples that are outside their
training data set that are not guaranteed to be accurate [36] and are unable
to identify these predictions as being OOD. As such, we incorporated the
desired qualities of the GP into the DNN by adding a fixed size lower rank
approximation of the GP with an RBF kernel, $K=\Phi\Phi^{T}$, at the final
layer using random fourier features (RFF) as defined in [37]. Although this
provides a uncertainty estimation that is distance aware, it is not distance
preserving because there is no guarantee that distance between the input data
is preserved at the hidden layer where the GP approximation is applied. In
order to make the DNN distance preserving, we used the bi-Lipschitz constraint
as part of the training loss function, as shown in Equation 6.

$L_{1}\times||x_{1}-x_{2}||\leq||h_{x_{1}}-h_{x_{2}}||\leq
L_{2}\times||x_{1}-x_{2}||$ (6)

Here $x$ is the input feature vector and $h_{x}$ be the last hidden layer
output. To summarize, we implemented the same model architecture used for the
other methods in this paper, however, we introduced a GP RBF kernel
approximation with 256 RFF and modified the loss function to ensure the
distance between input and hidden layers are preserved using soft bi-Lipschitz
constraint using $L_{1}=0.75$ and $L_{2}=1.25$.

##  IV RESULTS

Traditional deep learning models are deterministic and provide a prediction
for each input with no measure of confidence associated with the prediction.
Providing methods with reliable predictive uncertainty for ML models is
critical for real-world applications. As discussed in the previous section,
there are number of methods being proposed in the literature to make the DL
models uncertainty aware. In this section we compare the performance of the
methods presented in Section III for in-distribution and out-of-distribution
samples as it applies to the prediction for the FNAL Booster accelerator
complex. The input samples that are independent and identically distributed
(iid) as training data (in-distribution), if the underlying system/data
produces noisy labels, a DL model will learn to produce the mean. We expect
the uncertainty values for such predictions to reliably represent the variance
in the underlying data labels. The input samples that are dissimilar to the
training samples, called OOD samples, are most difficult for a DL model to
provide the prediction accurately; most of the time, the model will produce
inaccurate predictions leading to unreliable results. A prediction without the
associated calibrated, distance-aware uncertainty quantification results in a
system without contextual information required to select a safe response for a
prediction. We require the uncertainty values for OOD predictions to be high
indicating low confidence in the prediction.

To compare the results, we trained all the models with similar architecture
and data sets. The models consist of three convolutional blocks where each
block contained a 1-dimensional convolutional layer with 32 filters of size 3
followed by a maxpooling layer, and dropout layer with probability of 0.1. The
output of the third convolutional block is then flattened to process through a
dense layer containing 256 nodes leading to the output layer. The differences
between the three models are in how they quantify uncertainty. The DGPA model
has a Gaussian approximation layer as the output layer, BNN has a vanilla
dense layer as output layer but the dropouts are kept on during inference, and
DQR has multiple dense layers to produce output for different quantiles.

For this study, the raw data were processed using MinMax scaling and
restructured so that 15 previous timesteps of the input variables were fed
into the models to predict the next timestep forward in the output variables.
We divided the data into orthogonal samples, $80\%$ for training and $20\%$
for testing. The samples were further filtered by explicitly excluding
contributions when the main injector lower bound current (I:IB) had a value
that exceeded 0.995. This filter was used to create a in distribution only
training samples that would prevent the models to see the cyclic high
amplitude in the predicted variable (B:VIMIN). The relationship between the
filtered variable (I:IB) and the predicted variable (B:VIMIN) is shown in
Figure 3.

Figure 3: Comparison between the main injector lower bend current and the
compensated minimum GMPS current.

###  IV.1 In-distribution results

For in-distribution input samples, deep learning models are expected to have
accurate predictions with an uncertainty estimation that is consistent with
the training data. To compare the predictive performance along with the
uncertainty quantification for these models we use a set of standard metrics
including R-square, Root Mean Square Error (RMSE) between the ground truth
labels and the predictions, Mean Absolute Calibration Error (MACE), and Root
Mean Square Calibration Error (RMSCE) from the uncertainty toolkit [38]. Table
2 shows the values of these metrics for each of the three models before
performing any calibration.

All three models have very similar predictive performance in terms of $R^{2}$
and RMSE, however, their uncertainty estimations results vary. The DQR and
DGPA models produced similar uncertainty values that accurately represent the
standard deviation in the underlying training sample. The BNN model, on the
other hand, produced a MACE and RMSCE that is an order of magnitude higher
than the other two methods prior to offline calibration. This is also evident
in Figure 4, the BNN model underestimate the uncertainties as compared to DGPA
and DQR which produces accurate uncertainty estimates. The uncertainty
estimation from the BNN before calibration can vary significantly since the
dropout fraction during inference has yet to be optimized. However, this
illustrates the fact that this method cannot be used in real-time scenarios
where post-training calibration is not an option.

Table 2: In distribution prediction performance for BNN, DQR, and DGPA models.

Model | $R^{2}$ | RMSE | MACE | RMSCE
---|---|---|---|---
BNN | 0.825 | 0.033 | 0.224 | 0.251
DQR | 0.907 | 0.024 | 0.015 | 0.017
DGPA | 0.856 | 0.030 | 0.016 | 0.018

Figure 4: Comparison of miscalibration between DGPA, BNN and DQR for in-
distribution data. The shaded area represents the amount of miscalibration in
the uncertainty estimation with respect to the true labels.

###  IV.2 Out-of-distribution results

The majority of ML applications assumes that the test samples for a model are
iid and similar to the training data. Unfortunately, in practice this
assumption doesn’t always hold. When a test data draws from an out of training
distribution sample, the trained model is not guaranteed to produce accurate
predictions[36]. Providing an uncertainty estimation consistent with the
distance between in-distribution and out-of-distribution data is desirable. We
expect the model to produce high uncertainty for inaccurate predictions and
lower uncertainty when the predictions are accurate. We expect this
relationship between uncertainty and error even for the OOD samples. To
evaluate the ability of each model’s technique to estimate the OOD
uncertainty, we created two scenarios detailed below.

####  IV.2.1 Scenario 1

For the first scenario, we trained the models with the in-distribution samples
and evaluated their performance using the full test sample.

Figure 5: Side by side comparison of the predictive performance as well as
uncertainty quantification results for in-distribution and out-of-distribution
samples for DGPA, BNN, and DQR respectively. The middle region with the high
frequency component on the time series represent OOD samples while the initial
and tail-end regions represent in-distribution data samples.

The top three plots in Figure 5 show each model’s predictions and
uncertainties using data that includes the OOD samples. As can be seen, the
predictions from all three models degrade in the OOD region, which is
expected. The respective uncertainty values are expected to correlate to the
deviation from the ground truth labels. Table 3 describes the same set of
metrics used for in-distribution evaluation but for OOD samples. The $R^{2}$
and RMSE is similar for all three models. However, the uncertainty estimation
varies. Both BNN and DQR underestimate the OOD uncertainties as show in Figure
6.

In contrast, the DPGA has approximately 3x smaller MACE and RMSCE than the DQR
and BNN which indicated a better overall prediction uncertainty estimation.
These results are without any offline calibration.

Figure 6: Comparison of miscalibration between DGPA, BNN and DQR for OOD data.
The shaded area represents the amount of miscalibration in the uncertainty
estimation with respect to the true labels Table 3: Out of distribution
prediction performance for BNN, DQR, and DGPA models.

Model | $R^{2}$ | RMSE | MACE | RMSCE
---|---|---|---|---
BNN | 0.700 | 0.103 | 0.252 | 0.282
DQR | 0.836 | 0.076 | 0.275 | 0.309
DGPA | 0.784 | 0.088 | 0.091 | 0.100

####  IV.2.2 Scenario 2

Figure 7 shows the second scenario where we monotonically increase one of the
key input variable VIMIN until the data samples enters into a region of
feature space that is out of training distribution. VIMIN is chosen because it
is one of the key variable affecting our target. These data samples are then
fed into the models for inference and uncertainty quantification. Since these
data samples are not within the training distribution they are seen as OOD by
the models. Similar to the above discussed OOD results, the uncertainty values
from the BNN and DQR are underestimated as compared to DGPA when the data
samples goes into OOD region. From Figure 5 and Figure 6 as well as the Figure
7 it is clear that for OOD samples, DGPA produces more accurate uncertainty
estimates as compared to BNN and DQR, both of which underestimate the
uncertainty values.

Figure 7: Comparison of the predictive performance with uncertainty
quantification between DGPA, BNN, and DQR respectively for a manually induced
OOD. The vertical line at 450 indicates the area where the VIMIN input
variable was incrementally increasing to induce the OOD scenario.

##  V Summary and Conclusions

In this paper, we compared three different DNN techniques that estimate the
prediction uncertainty. We present the DGPA technique as a new approach to
estimating prediction uncertainties; DGPA is self-calibrated and includes an
awareness of out-of-distribution uncertainty. The results show that all models
provide similar performance for the predictions for in-distribution and out-
of-distribution from Scenario 1. The in-distribution uncertainty estimation
for the DQR and DGPA models provides excellent results and is significantly
better than the BNN model before calibration. Additionally, the DGPA provides
the best uncertainty estimation for the out-of-distribution predictions
described in Scenario 1. Although we cannot quantify the explicit expected
uncertainty for Scenario 2, we can see that the uncertainty estimation from
the DPGA is larger than the other two methods where the DQR provides a very
small uncertainty estimation. In conclusion, the results from this study on
the FNAL Booster Accelerator Complex data suggestion that the DGPA model
provides the best single inference calibrated model for in and out of
distribution uncertainty estimation for all scenarios. This was achieved by
using a fixed size GP RBF kernel approximation and applying the bi-Lipschitz
constraint in the loss function. Additional research should be conducted to
better understand the trade-off from kernel approximation size and how this
can be applied to real-time systems where the hardware could be a constraint.

###### Acknowledgements.

This research was supported by the U.S. Department of Energy, through the
Office of Advanced Scientific Computing Research’s “Data-Driven Decision
Control for Complex Systems (DnC2S)” project. Pacific Northwest National
Laboratory is operated by Battelle Memorial Institute for the U.S. Department
of Energy under Contract No. DE-AC05-76RL01830. Oak Ridge National Laboratory
is operated by UT-Battelle LLC for the U.S. Department of Energy under
contract number DE-AC05-00OR22725. The Jefferson Science Associates (JSA)
operates the Thomas Jefferson National Accelerator Facility for the U.S.
Department of Energy under Contract No. DE-AC05-06OR23177.

## References

  * Amundson  _et al._ [2014] J. Amundson, A. Macridin, and P. Spentzouris, High performance computing modeling advances accelerator science for high energy physics, IEEE Comput. Sci. Eng. 16, 32 (2014).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, Nature 521, 436 (2015).
  * Goodfellow  _et al._ [2016] I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (MIT Press, 2016) http://www.deeplearningbook.org.
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Rev. Mod. Phys. 91, 045002 (2019), arXiv:1903.10563 .
  * DeepMind [2018] DeepMind, Safety-first AI for autonomous data centre cooling and industrial control, https://deepmind.com/blog/safety-first-ai-autonomous-data-centre-cooling-and-industrial-control/ (2018).
  * Davies [2018] A. Davies, The WIRED Guide to Self-Driving Cars, https://www.wired.com/story/guide-self-driving-cars/ (2018).
  * Jennifer Langston [2018] Jennifer Langston, How AI is building better gas stations and transforming Shell’s global energy business, https://blogs.microsoft.com/ai/shell-iot-ai-safety-intelligent-tools/ (2018).
  * Will Knight [2018] Will Knight, This Factory Robot Learns a New Job Overnight, https://www.technologyreview.com/s/601045/this-factory-robot-learns-a-new-job-overnight/ (2018).
  * Emma  _et al._ [2018] C. Emma, A. Edelen, M. J. Hogan, B. O’Shea, G. White, and V. Yakimenko, Machine learning-based longitudinal phase space prediction of particle accelerators, Phys. Rev. Accel. Beams 21, 112802 (2018).
  * Sanchez-Gonzalez  _et al._ [2017] A. Sanchez-Gonzalez _et al._ , Machine learning applied to single-shot x-ray diagnostics in an XFEL, Nature Commun. 8, 15461 (2017), arXiv:1610.03378 [physics.data-an] .
  * Wielgosz  _et al._ [2017] M. Wielgosz, A. Skoczeń, and M. Mertik, Using LSTM recurrent neural networks for monitoring the LHC superconducting magnets, Nucl. Instrum. Meth. A 867, 40 (2017), arXiv:1611.06241 [physics.ins-det] .
  * Scheinker and Gessner [2015] A. Scheinker and S. Gessner, Adaptive method for electron bunch profile prediction, Phys. Rev. ST Accel. Beams 18, 102801 (2015).
  * Scheinker  _et al._ [2020] A. Scheinker, C. Emma, A. L. Edelen, and S. Gessner,  _Advanced Control Methods for Particle Accelerators (ACM4PA) 2019 Workshop Report_, Tech. Rep. (2020) arXiv:2001.05461 .
  * Scheinker  _et al._ [2018] A. Scheinker, A. Edelen, D. Bohler, C. Emma, and A. Lutman, Demonstration of model-independent control of the longitudinal phase space of electron beams in the linac-coherent light source with femtosecond resolution, Phys. Rev. Lett. 121, 044801 (2018).
  * Blokland  _et al._ [2021] W. Blokland, P. Ramuhalli, C. Peters, Y. Yucesan, A. Zhukov, M. Schram, K. Rajput, and T. Jeske, Uncertainty aware anomaly detection to predict errant beam pulses in the SNS accelerator, (2021), arXiv:2110.12006 [physics.acc-ph] .
  * Rescic  _et al._ [2020] M. Rescic, R. Seviour, and W. Blokland, Predicting particle accelerator failures using binary classifiers, Nucl. Instrum. Meth. A 955, 163240 (2020).
  * Reščič  _et al._ [2022] M. Reščič, R. Seviour, and W. Blokland, Improvements of pre-emptive identification of particle accelerator failures using binary classifiers and dimensionality reduction, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 1025, 166064 (2022).
  * Tennant  _et al._ [2020] C. Tennant, A. Carpenter, T. Powers, A. Shabalina Solopova, L. Vidyaratne, and K. Iftekharuddin, Superconducting radio-frequency cavity fault classification using machine learning at jefferson laboratory, Phys. Rev. Accel. Beams 23, 114601 (2020).
  * Powers and Solopova [2019] T. Powers and A. Solopova, CEBAF C100 Fault Classification based on Time Domain RF Signals, in  _19th International Conference on RF Superconductivity (SRF 2019)_ (2019) p. WETEB3.
  * St. John  _et al._ [2021] J. St. John, C. Herwig, D. Kafkes, J. Mitrevski, W. A. Pellico, G. N. Perdue, A. Quintero-Parra, B. A. Schupbach, K. Seiya, N. Tran, M. Schram, J. M. Duarte, Y. Huang, and R. Keller, Real-time artificial intelligence for accelerator control: A study at the fermilab booster, Phys. Rev. Accel. Beams 24, 104601 (2021).
  * Kafkes and Schram [2021] D. Kafkes and M. Schram, Developing Robust Digital Twins and Reinforcement Learning for Accelerator Control Systems at the Fermilab Booster, in  _12th International Particle Accelerator Conference_ (2021) arXiv:2105.12847 [physics.acc-ph] .
  * Edelen  _et al._ [2016] A. L. Edelen, S. G. Biedron, B. E. Chase, D. Edstrom, S. V. Milton, and P. Stabile, Neural Networks for Modeling and Control of Particle Accelerators, IEEE Trans. Nucl. Sci. 63, 878 (2016), arXiv:1610.06151 [physics.acc-ph] .
  * Hirlaender and Bruchon [2020] S. Hirlaender and N. Bruchon, Model-free and bayesian ensembling model-based deep reinforcement learning for particle accelerator control demonstrated on the fermi fel (2020), arXiv:2012.09737 .
  * Mishra  _et al._ [2021] A. A. Mishra, A. Edelen, A. Hanuka, and C. Mayes, Uncertainty quantification for deep learning in particle accelerator applications, Phys. Rev. Accel. Beams 24, 114601 (2021).
  * Minorsky. [1922] N. Minorsky., Directional stability of automatically steered bodies, J. Am. Soc. Nav. Engineers 34, 280 (1922).
  * Ziegler and Nichols [1942] J. G. Ziegler and N. B. Nichols, Optimum settings for automatic controllers, J. Dyn. Syst. Meas. Control 64, 759 (1942).
  * Kafkes and St. John [2021] D. Kafkes and J. St. John, BOOSTR: A dataset for accelerator control systems, Data 6, 42 (2021), arXiv:2101.08359 .
  * Gal and Ghahramani [2016] Y. Gal and Z. Ghahramani, Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in  _international conference on machine learning_ (PMLR, 2016) pp. 1050–1059.
  * Koenker [2005] R. Koenker,  _Quantile Regression_, Econometric Society Monographs (Cambridge University Press, 2005).
  * Rahimi and Recht [2007] A. Rahimi and B. Recht, Random features for large-scale kernel machines, in  _Advances in Neural Information Processing Systems_, Vol. 20, edited by J. Platt, D. Koller, Y. Singer, and S. Roweis (Curran Associates, Inc., 2007).
  * Rasmussen [2004] C. E. Rasmussen, Gaussian processes in machine learning, in  _Advanced Lectures on Machine Learning: ML Summer Schools 2003, Canberra, Australia, February 2 - 14, 2003, Tübingen, Germany, August 4 - 16, 2003, Revised Lectures_, edited by O. Bousquet, U. von Luxburg, and G. Rätsch (Springer Berlin Heidelberg, Berlin, Heidelberg, 2004) pp. 63–71.
  * Liu  _et al._ [2020a] J. Liu, Z. Lin, S. Padhy, D. Tran, T. Bedrax Weiss, and B. Lakshminarayanan, Simple and principled uncertainty estimation with deterministic deep learning via distance awareness, in  _Advances in Neural Information Processing Systems_, Vol. 33, edited by H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin (Curran Associates, Inc., 2020) pp. 7498–7512.
  * Wilson  _et al._ [2015] A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing, Deep kernel learning (2015).
  * Liu  _et al._ [2020b] J. Z. Liu, Z. Lin, S. Padhy, D. Tran, T. Bedrax Weiss, and B. Lakshminarayanan, Simple and principled uncertainty estimation with deterministic deep learning via distance awareness (2020b).
  * Abdar  _et al._ [2021] M. Abdar, F. Pourpanah, S. Hussain, D. Rezazadegan, L. Liu, M. Ghavamzadeh, P. Fieguth, X. Cao, A. Khosravi, U. R. Acharya, _et al._ , A review of uncertainty quantification in deep learning: Techniques, applications and challenges, Information Fusion 76, 243 (2021).
  * Cortes  _et al._ [2009] C. Cortes, M. Mohri, and A. Rostamizadeh, Learning non-linear combinations of kernels, in  _Advances in Neural Information Processing Systems_, Vol. 22, edited by Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta (Curran Associates, Inc., 2009).
  * Liu  _et al._ [2020c] J. Z. Liu, Z. Lin, S. Padhy, D. Tran, T. Bedrax-Weiss, and B. Lakshminarayanan, Simple and principled uncertainty estimation with deterministic deep learning via distance awareness (2020c).
  * Chung  _et al._ [2021] Y. Chung, I. Char, H. Guo, J. Schneider, and W. Neiswanger, Uncertainty toolbox: an open-source library for assessing, visualizing, and improving uncertainty quantification, arXiv preprint arXiv:2109.10254 (2021).


# A Stackelberg viral marketing design for two competing players∗††thanks:
∗This work was supported by ANR through the grant NICETWEET No.
ANR-20-CE48-0009 and the PNRR project DECIDE No. 760069.

O. Lindamulage De Silva1, V. S. Varma1, M. Cao2, I.-C. Morărescu1,3, S.
Lasaulce 1
1 Université de Lorraine, CNRS, CRAN, F-54000 Nancy, France,
<EMAIL_ADDRESS>Faculty of Mathematics and
Natural Sciences, ITM, University of Groningen, The Netherlands3associated
with Automation Department, Technical University of Cluj-Napoca,
Memorandumului 28, 400114 Cluj-Napoca, Romania

###### Abstract

A Stackelberg duopoly model in which two firms compete to maximize their
market share is considered. The firms offer a service/product to customers
that are spread over several geographical regions (e.g., countries, provinces,
or states). Each region has its own characteristics (spreading and recovery
rates) of each service propagation. We consider that the spreading rate can be
controlled by each firm and is subject to some investment that the firm does
in each region. One of the main objectives of this work is to characterize the
advertising budget allocation strategy for each firm across regions to
maximize its market share when competing. To achieve this goal we propose a
Stackelberg game model that is relatively simple while capturing the main
effects of the competition for market share. By characterizing the strong/weak
Stackelberg equilibria of the game, we provide the associated budget
allocation strategy. In this setting, it is established under which conditions
the solution of the game is the so-called “winner takes all”. Numerical
results expand upon our theoretical findings and we provide the equilibrium
characterization for an example.

###### Index Terms:

Winner takes all, viral marketing, resource allocation, Stackelberg solution.

##  I Introduction

Viral marketing (VM) is a strategy to promote different products/services
using social networks. The information spreads as a virus from one person to
their family, friends and colleagues. To model the spread of information or
services within a population, various mathematical frameworks have been
developed. In the context of viral marketing, two prominent modeling
approaches are opinion dynamics models (see [1, 2] for instance) and epidemic
models (see [3, 4] for instance). Opinion dynamics models focus on capturing
how opinions and preferences evolve among individuals in a social network,
considering factors such as interaction pattens, biased influences, and
coupled decision-making processes. On the other hand, epidemic models capture
the transmission of diseases or information within a population, incorporating
parameters such as infection rates, recovery rates, and population
connectivity, which play a crucial role in understanding the dynamics of
service propagation and devising effective marketing strategies to maximize
the spread and adoption of the service. In this paper, we are not focusing on
a specific dynamic but assume that the equilibrium in each region corresponds
to the so-called “winner takes all” [5], for which a particular case of such
model is the steady-state of a multi-virus SIS system [6]. It is noteworthy
that such an equilibrium has been observed in real life, see for instance the
classical example of Facebook and Myspace [6]. In the context of viral
marketing duopoly, the concept of “winner takes all” refers to a scenario
where two competing firms engage in viral marketing campaigns to gain a larger
market share. In this situation, the firm that effectively leverages viral
marketing techniques and achieves widespread adoption of its product or
service tends to capture the majority of the market, leaving the competitor
with a significantly smaller share.

It has been proven that targeted marketing combined with social network
spreading has advantages over conventional mass campaigns, including cost-
effectiveness and the ability to reach specific customer groups [2].
Basically, the authors formulated the problem as an optimal budget allocation
and they shown that most individuals have to be targeted by the marketing
campaign in order to get a maximum profit. The problem of competition to get a
larger market share has been addressed in [7]. The authors introduced a
duopoly model which accounts for the knowledge of opinion dynamics through a
social network and characterized the Nash strategies of the players. Unlike
these works that take advantage of the node centrality (network topology) but
rely on linear opinion dynamic models for the service spreading, we assume
here a winner takes all model based on [6]. This setup is more suitable for
certain types of products/services, such as video streaming and activities in
social networking platforms.

The main goal of this paper is to formulate a two-level Stackelberg game and
characterize its solutions. The analysis does not rely on a specific dynamic
for the information spreading but assumes that at the equilibrium the “winner
takes all”. We point out that the existing literature on multi-virus epidemic
emphasize the aforementioned type of equilibrium [6, 3, 4]. We consider that
two players compete over several regions to get a higher market share. The
budget allocated by a firm to a certain region modifies the spreading rate of
the associated service in that region. Note that we are analyzing the steady
state revenue of each firm. This combined with the “winner takes all” behavior
allows us to decouple the analysis of the investment strategy for each region.
At the same time, we highlight that the design of the budget allocation
strategy is subject to overall fixed budget constraints.

Note that each player has to solve a budget allocation problem which is
different from the ones that can be found in the literature. In [8], a dynamic
optimization problem under budget constraints is formulated to control a
single-virus SIS model. In [9], optimal control of joint multi-virus infection
and information dissemination is considered for the sensitive-warned-infected-
recovered-susceptible (SWIRS) model, without budget constraints. The
differences between the present work and the existing results on epidemic
control (e.g., [10, 11, 12, 13, 14]) are mainly related to the fact that our
model can handle a multi-virus SIS epidemic model by considering budget
constraints. Two existing works are relatively close to the present one. The
first is [15] in which the authors formulate a static and strategic form game
to deal with a bi-virus SIS epidemic model over a single region without a
budget constraint. The second one is [16] in which only one player solves the
optimal budget allocation problem over several regions.

This paper is structured as follows. Sec. II provides the problem formulation,
going from the VM model up to the Stackelberg game analyzed in the paper. The
main result of this work is presented at the end of Sec. II. Sec. III is
devoted to the proof of the main result. A numerical example illustrates our
theoretical findings in Sec. IV and some concluding remarks are provided in
Sec. V.

##  II Problem statement

We consider a set of two firms competing over $K>1$ regions (e.g., countries,
provinces, or states) to maximize their market share. Let
$\mathcal{K}:=\\{1,\ldots,K\\}$ and $\mathcal{M}:=\\{1,2\\}$ be the set of
regions and firms, respectively. For a given region $k\in\mathcal{K}$ and a
firm $m\in\mathcal{M}$, we respectively denote by (i) $\gamma_{mk}$ the
spreading rate of the service of Firm $m$ in the region $k$; (ii)
$\delta_{mk}$ the churn rate at which the individuals from Region $k$ decide
to dispense with the services of Firm $m$. In this work, we assume that there
is a simple linear relationship between allocated budget and $\gamma_{mk}$
dissemination rate, but more advanced and faithful models could be
investigated in future work. In practice, the spread rate is influenced by
numerous factors and it would be more accurate to assume that firms focus on
captivating content and identifying influential users or communities to
indirectly influence the spread rate of viral marketing campaigns [17, 18,
19]. On top of this, we assume that Firm $m$ has a given advertising budget
$B_{m}$ to allocate between the $K$ regions in order to maximize the number of
its subscribers. We denote by
$\gamma_{m}:=(\gamma_{m,1},\ldots,\gamma_{mk})\in\mathbb{R}_{\geq 0}^{K}$ the
action vector of Firm $m$ such that $\sum_{k=1}^{K}\gamma_{mk}\leq B_{m}$ and
$\gamma:=(\gamma_{1},\gamma_{2})\in\mathbb{R}_{\geq 0}^{2K}$ the action
profile. We also make a slight notation abuse by using the following notation:
$\gamma:=(\gamma_{m},\gamma_{-m})$ for $m\in\mathcal{M}$.

Firm $1$ will be referred to as the _leader_ and Firm $2$ as the _follower_.
This is because the leader acts on the network first and the follower acts
after a sufficiently long time such that it can observe and react to the
action made by the leader.

###  II-A Viral marketing model

By considering that the control action $\gamma_{mk}$ is a constant on the
working phase $[0,+\infty[$. Thus, the follower can be said to start
influencing nodes at time $0$. The fraction of individuals in Region
$k\in\mathcal{K}$ who subscribe to the services of Firm $m\in\mathcal{M}$ is
denoted by $x_{mk}\in[0,1]$. For the interest of practicality, we consider the
following assumption.

###### Assumption 1

Each individual subscribes to at most one service.

This assumption is justified by the consideration of average individual
behavior, which aligns with the notion that individuals typically opt for a
single subscription. The fraction of individuals in Region $k\in\mathcal{K}$
who have not subscribed to any services is denoted by $s_{k}\in[0,1]$. In what
follows we denote by $x_{mk}^{\infty}:=\lim\limits_{t\to\infty}x_{mk}(t)$
$\forall k\in\mathcal{K},\forall m\in\mathcal{M}$ and we suppose that:

$(x_{1k}^{\infty},x_{2k}^{\infty})=\left\\{\begin{array}[]{ll}\displaystyle\Big{(}1-\frac{\delta_{1k}}{\gamma_{1k}},0\Big{)}&\text{if
}\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}\geq\max\left(\displaystyle\frac{\gamma_{2k}}{\delta_{2k}},1\right)\\\
\displaystyle\Big{(}0,1-\frac{\delta_{2k}}{\gamma_{2k}}\Big{)}&\text{if
}\displaystyle\frac{\gamma_{2k}}{\delta_{2k}}\geq\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\\\
(0,0)&\text{otherwise.}\end{array}\right.$ (1)

The presence of $\pi>0$ in (1) is motivated by practical reasons well-known in
the economics literature. The leader has already established its market by the
time when the follower enters the market. Thus, the follower has to invest a
little more to convince the customers to switch services. In the economics
literature $\pi$ is called a barrier to entry. A barrier to entry is a
condition that makes it difficult for new firms to enter a market and compete
with established firms. Barriers can take various forms, such as economies of
scale, brand loyalty, access to distribution channels, patents and copyrights,
government regulations, and high capital requirements [20, 21]. The analysis
in [6, Section 4.2] and [22] provide a bi-virus (Susceptible-Infected-
Susceptible) SIS model applied for viral marketing, such that the system
converges to a point closed to (1). The main difference is that the case when
$\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}=\displaystyle\frac{\gamma_{2k}}{\delta_{2k}}$
leads to multiple equilibria, and is avoided in our model due to the barrier
to entry $\pi$. Additionally, by setting the churn rates $\delta_{mk}$ to
zero, the model corresponds to a Colonel Blotto game, which is commonly
employed in viral marketing literature [23].

###  II-B Game model

As previously stated, each Firm $m\in\mathcal{M}$ solves a budget allocation
problem that maximizes a revenue under the global budget (denoted by
${B}_{m}$) constraint i.e., a feasible strategy for Firm $m$ belong to the set
$\Gamma_{m}:=\\{\gamma_{m}\in\mathbb{R}_{\geq 0}^{K}:\
\sum_{k=1}^{K}\gamma_{mk}\leq\mathrm{B_{m}}\\}$. The budget ${B}_{m}$ is
imposed for practical reasons (e.g., a firm has a given finite investment
budget).

Similar to [24], we consider the utility of each Firm $m\in\mathcal{M}$, such
as given by

$\begin{array}[]{l}\displaystyle
u_{1}(\gamma_{1},\gamma_{2}):=\sum_{k=1}^{K}p_{1k}\left(1-\frac{\delta_{1k}}{\gamma_{1k}}\right)\leavevmode\resizebox{147.4292pt}{}{$\mathbb{1}_{\mathcal{H}_{1}(k,\gamma_{2k})}(\gamma_{1k})$}\\\
\displaystyle
u_{2,\pi}(\gamma_{2},\gamma_{1}):=\sum_{k=1}^{K}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)\leavevmode\resizebox{147.4292pt}{}{$\mathbb{1}_{\mathcal{H}_{2,\pi}(k,\gamma_{1k})}(\gamma_{2k})$}\end{array}$

(2)

where: (i) $\pi>0$ is the barrier to entry for the follower, which is fixed;
(ii) $p_{mk}$ the contribution of Region $k$ to the revenue of firm $m$ when
all its individuals subscribe to the service $m$; (iii) the quantity
$\big{(}1-\frac{\delta_{mk}}{\gamma_{mk}}\big{)}$ represents
$x_{mk}^{\infty}(\gamma,\pi)$ in view of (1) and $\mathbb{1}$ is the indicator
function, where

$\mathcal{H}_{1}(k,\gamma_{2k}):=\left\\{\gamma_{1k}\in\mathbb{R}_{\geq 0}:\
\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}\geq\max\left(\displaystyle\frac{\gamma_{2k}}{\delta_{2k}},1\right)\right\\},$

and

$\mathcal{H}_{2,\pi}(k,\gamma_{1k}):=\left\\{\gamma_{2k}\in\mathbb{R}_{\geq
0}:\
\displaystyle\frac{\gamma_{2k}}{\delta_{2k}}\geq\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\right\\}.$

In this paper, we analyze a static Stackelberg game with a leader (Firm 1) and
a follower (Firm 2). The utility of Firm 1 depends only on $\gamma_{1}$ as
Firm 2 reacts with a best response strategy. We adopt a pessimistic approach
for the leader, focusing on the analysis of weak Stackelberg equilibrium [25].
The leader’s utility function, considering a given barrier to entry $\pi>0$
for the follower, is formulated as:

$\left.\begin{array}[]{ll}\displaystyle
u_{1,\pi}^{S}(\gamma_{1}):=\min\limits_{\gamma_{2}\in\mathrm{BR}_{2,\pi}(\gamma_{1})}u_{1}(\gamma_{1},\gamma_{2}),\\\
\displaystyle\mathrm{BR}_{2,\pi}(\gamma_{1})=\operatorname*{\text{argmax}\,}_{\gamma_{2}\in\Gamma_{2}}u_{2,\pi}(\gamma_{2},\gamma_{1}).\end{array}\right.$
(3)

It should be noted that the analysis presented in this paper also applies when
considering an optimistic formulation for the leader, namely, strong
Stackelberg equilibrium with
$u_{1,\pi}^{S}(\gamma_{1})=\max\limits_{\gamma_{2}\in\mathrm{BR}_{2,\pi}(\gamma_{1})}u_{1}(\gamma_{1},\gamma_{2})$.
The goal of this work is to analyse a regular Stackelberg solution (S), for a
given $\pi>0$ of the game
$\mathcal{G}_{\pi}^{S}:=\left(\mathcal{M},(\Gamma_{m})_{1\leq m\leq
2},(u_{m,\pi}^{S})_{1\leq m\leq 2}\right)$ in which: the players are Firm $1$
and Firm $2$; the action space of Firm $m$ is given by $\Gamma_{m}$; the
individual utility function of each firm is given by $u_{1,\pi}^{S}$ in (3)
and $u_{2,\pi}^{S}=u_{2,\pi}$ in (2). Firm $m\in\mathcal{M}$ expresses its
interests by setting the potential revenues
$(p_{m,1},\ldots,p_{mk})\in\mathbb{R}_{\geq 0}^{K}$, whereas the set of action
$\Gamma_{m}$ is imposed by the capacity of investment of Firm $m$ given by the
budget $\mathrm{B_{m}}$. In addition, we emphasize that the theoretical
results established in this paper hold for a multistage game setup in which
the one-shot game is repeated at each stage (for which the parameters are
updated) and different constant control actions are applied during it.

Let us recall the definition of a weak Stackelberg solution. The strategy
${\gamma}_{\pi}^{\mathrm{S}}$ is a weak Stackelberg solution of the game
$\mathcal{G}_{\pi}^{S}$ if it is a solution of the system of equations:

$\left.\begin{array}[]{ll}\gamma_{1,\pi}^{\mathrm{S}}\in\operatorname*{\text{argmax}\,}\limits_{\gamma_{1}\in\Gamma_{1}}u_{1,\pi}^{S}(\gamma_{1})\text{
and
}\gamma_{2,\pi}^{\mathrm{S}}\in\operatorname*{\text{argmax}\,}\limits_{\gamma_{2}\in\Gamma_{2}}u_{2,\pi}^{S}(\gamma_{2},\gamma_{1,\pi}^{\mathrm{S}}).\end{array}\right.$

###  II-C Main result

As mentioned above, we investigate the weak Stackelberg solution of the game
$\mathcal{G}_{\pi}^{S}$. In what follows, we denote by:

$\mathrm{BR}_{2,\pi}^{-}(\gamma_{1}):=\operatorname*{\text{argmin}\,}\limits_{\gamma_{2}\in\mathrm{BR}_{2,\pi}(\gamma_{1})}u_{1}(\gamma_{1},\gamma_{2})$

the best response of Firm $2$ which minimizes the profit of the Firm $1$ when
it uses the strategy $\gamma_{1}\in\Gamma_{1}$. Let
$\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1}):=\mathrm{P}_{k}(\mathrm{BR}_{2,\pi}^{-}(\gamma_{1}))$
where $\forall\gamma_{2}\in\Gamma_{2},\ k\in\mathcal{K}$ one has
$\mathrm{P}_{k}(\gamma_{2})=\gamma_{2k}$.

Let $\mathcal{K}_{1}\in 2^{K}$ representing a given set of regions of
investment under consideration by Firm 1, let us define the following set

$\begin{array}[]{ll}\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1}):=\Bigg{\\{}\gamma_{1}\in\mathbb{R}_{\geq
0}^{K}\Big{|}\\\ \displaystyle\forall k\in\mathcal{K}_{1},\
\frac{\gamma_{1k}}{\delta_{1k}}\geq\max\left(\frac{\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1})}{\delta_{2k}},1\right),\\\
\displaystyle\forall k\in\mathcal{K}\setminus\mathcal{K}_{1},\
\frac{\gamma_{1k}}{\delta_{1k}}\leq\max\left(\frac{\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1})}{\delta_{2k}}-\pi,1\right)\Bigg{\\}},\end{array}$

and $\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1}):=$

$\left\\{\begin{array}[]{ll}&\sum_{k\in\mathcal{K}_{1}}p_{1k}\left(\displaystyle
1-\frac{\delta_{1k}}{\gamma_{1k}}\right)\text{ if
}\gamma_{1}\in\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})\neq\emptyset,\\\
&-1\hskip 88.2037pt\text{ otherwise. }\end{array}\right.$ (4)

In simpler terms, the expression $\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})$
represents the set of investment options available to Firm 1 when Firm 2
applies a best response strategy. This set allows Firm 1 to determine the
regions and values of investment that it can choose from. On the other hand,
$\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})$ represents the corresponding
revenue that Firm 1 can generate based on its chosen investment strategy
$\mathcal{K}_{1}$ and spread rate $\gamma_{1}$. We are now ready to state the
main result of this paper, which will be proven in the following section.

###### Theorem 1

The Stackelberg strategy $\gamma_{\pi}^{S}$ is obtained by solving:

$(\mathcal{K}_{1,\pi}^{S},\gamma_{1,\pi}^{S}):=\operatorname*{\text{argmax}\,}_{\mathcal{K}_{1}\in
2^{\mathcal{K}}}\operatorname*{\text{argmax}\,}_{\gamma_{1}\in\Gamma_{1}}\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})$
(5)

where $\widehat{\Gamma}_{1,\pi}(\mathcal{K})$ is proven to be a convex set
(see proof) for any $\mathcal{K}$ and the maximization problem
$\max\limits_{\gamma_{1}\in\Gamma_{1}}\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})$
is a convex OP for a given $\mathcal{K}_{1}$. The follower’s strategy at the
Stackelberg equilibria is given by
$\gamma_{2k,\pi}^{S}=\frac{\displaystyle\sqrt{p_{2k}\delta_{2k}}\displaystyle}{\displaystyle\sum_{\ell\in{\mathcal{K}}_{2,\pi}^{S}}\sqrt{p_{2\ell}\delta_{2\ell}}}\mathrm{B_{2}}$
if $k\in{\mathcal{K}}_{2,\pi}^{S}$ and $\gamma_{2k,\pi}^{S}=0$ otherwise,
where

$\begin{array}[]{l}\displaystyle\mathcal{K}_{2,\pi}^{S}:=\operatorname*{\text{argmax}\,}_{\mathcal{K}_{2}\in
2^{\mathcal{K}}}\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\displaystyle\sqrt{\frac{\delta_{2k}}{p_{2k}}}\sum_{\ell\in\mathcal{K}_{2}}\sqrt{p_{2\ell}\delta_{2\ell}}}{\mathrm{B_{2}}}\right)\\\
\hskip 71.13188pt\text{s.t}.\
\mathcal{K}_{2}\cap\mathcal{K}_{1,\pi}^{S}=\emptyset.\end{array}$

###### Remark 1

In the perspective of studying the strong Stackelberg solution, it is enough
to consider that the follower applies the strategy in
$\mathrm{BR}_{2,\pi}^{+}(\gamma_{1})$ where

$\mathrm{BR}_{2,\pi}^{+}(\gamma_{1}):=\operatorname*{\text{argmax}\,}\limits_{\gamma_{2}\in\mathrm{BR}_{2,\pi}(\gamma_{1})}u_{1}(\gamma_{1},\gamma_{2}).$

###### Remark 2

Although the optimization problem specified in Equation (5) theoretically
involves an exponential search space of $2^{K}$ regions. On one hand we
provide just a methodology to solve such problems (that does not involve
online solving and therefore computation time is not very important) and on
the other hand, in many practical scenarios, the number of regions $K$ in
which firms compete is relatively small, typically around $10$. Basically we
can consider that USA, Europe, South America represent regions for competition
at global level. If the competition is specific to a country than the number
of regions will be also relatively small. This manageable number of regions
significantly reduces the computational complexity associated with solving the
optimization problem. A key result introduced in Theorem 1 is by ensuring that
the maximization of $\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})$, for a
given $\mathcal{K}_{1}$ and for $\gamma_{1}\in\Gamma_{1}$ can be effectively
resolved through the utilization of a convex optimization algorithm.

In order to prove the main result of this paper, we first characterize the
best response strategy of Firm $2$ and finally characterize the weak
Stackelberg solution for a fixed barrier to entry $\pi>0$.

##  III Stackelberg strategy design

In this section, we propose a convex reformulation of the problem in (3). This
allows us to obtain the optimal budget allocation solution for the follower
and thus characterize the weak Stackelberg solution of the game
$\mathcal{G}_{\pi}^{S}$.

For a given leader’s strategy $\gamma_{1}\in\Gamma_{1}$ we denote by
$\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})$ as the set of
strategy of Firm $2$ to win the marketing battles in regions of
$\mathcal{K}_{2}\subseteq\mathcal{K}$ i.e.,
$\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1}):=\Bigg{\\{}\gamma_{2}\in\mathbb{R}_{\geq
0}^{K}\Bigg{|}$

$\begin{array}[]{ll}\displaystyle\forall k\in\mathcal{K}_{2},\
\frac{\gamma_{2k}}{\delta_{2k}}\geq\max\left(\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right),\\\
\displaystyle\forall k\in\mathcal{K}\setminus\mathcal{K}_{2},\
\frac{\gamma_{2k}}{\delta_{2k}}\leq\max\left(\frac{\gamma_{1k}}{\delta_{1k}},1\right)\Bigg{\\}}.\end{array}$
(6)

###  III-A Follower’s OP reformulation

The goal of this section is to reformulate the optimisation problem for the
follower ($m=2$) introduced in (3).

###### Proposition 1

Let $\gamma_{1}\in\Gamma_{1}$. The initial OP in (3) can be reformulated such
as
$\displaystyle\max_{\gamma_{2}\in\Gamma_{2}}u_{2,\pi}^{S}(\gamma_{2},\gamma_{1})=$

$\displaystyle\max_{\mathcal{K}_{2}\in
2^{\mathcal{K}}}\left[\begin{array}[]{ll}&\displaystyle\max_{\gamma_{2}}\displaystyle\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)\\\
&\hskip
28.45274pt\text{s.t.}\displaystyle\sum_{k=1}^{K}\gamma_{2k}\leq\mathrm{B_{2}}\\\
&\hskip
42.67912pt\gamma_{2}\in\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1}).\end{array}\right]$
$\hskip 48.36958pt\text{s.t.
}\exists\gamma_{2}\in\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1}):\
\sum_{k=1}^{K}\gamma_{2k}\leq\mathrm{B_{2}}$

###### Remark 3

It appears that, the condition
$``\exists\gamma_{2}\in\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1}):\
\sum_{k=1}^{K}\gamma_{2k}\leq\mathrm{B_{2}}^{\prime\prime}$ is verified, if
and only if
$\displaystyle\sum_{k\in\mathcal{K}_{2}}\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\leq\mathrm{B_{2}}$
i.e., the less restrictive strategy of
$\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})$ verifies the budget
constraint.

_$\mathbf{Proof.}$_ In view of (3), the follower’s OP can be rewritten in the
following manner by considering both cases of Firm $2$ winning or not in
region $k$:

$\displaystyle\max_{\gamma_{2}\in\Gamma_{2}}\sum_{\mathcal{K}_{2}\in
2^{\mathcal{K}}}\left[\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)\right]\mathbb{1}_{\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})}(\gamma_{2}),$

where $\widetilde{\Gamma}_{2,\pi}$ is introduced in (6) and $\mathbb{1}$ is
the indicator function.

Let $\gamma_{2},\widehat{\gamma}_{2}\in\Gamma_{2}$,
$\mathcal{K}_{2},\widehat{\mathcal{K}}_{2}\in 2^{\mathcal{K}}$ such that
$\mathcal{K}_{2}\neq\widehat{\mathcal{K}}_{2}$ and
$\mathbb{1}_{\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})}(\gamma_{2})=\mathbb{1}_{\widetilde{\Gamma}_{2,\pi}(\widehat{\mathcal{K}}_{2},\gamma_{1})}(\widehat{\gamma}_{2})=1$.
Then,
$\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})\cap\widetilde{\Gamma}_{2,\pi}(\widehat{\mathcal{K}}_{2},\gamma_{1})=\emptyset$
since $\mathcal{K}_{2}\neq\widehat{\mathcal{K}}_{2}$ and (III-A) can be
reformulated such as

$\max\limits_{\mathcal{K}_{2}\in
2^{\mathcal{K}}}\Bigg{[}\displaystyle\max_{\gamma_{2}\in\Gamma_{2}}\Big{[}\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)\Big{]}\mathbb{1}_{\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})}(\gamma_{2})\Bigg{]}.$

Furthermore, for a given $\mathcal{K}_{2}$ there exits a strategy
$\gamma_{2}\in\widetilde{\Gamma}_{m}(\mathcal{K}_{2},\gamma_{1})$ that
verifies the budget constraint if and only if the set $\mathcal{K}_{2}$
verifies
$\displaystyle\sum_{k\in\mathcal{K}_{2}}\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\leq\mathrm{B_{2}}$
i.e., the less restrictive action of
$\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1})$ verifies the budget
constraint. Hence, by using the indicator function and by adding the new
constraint on the feasibility sets of $\mathcal{K}_{2}$, the problem (III-A)
is equivalent to (1). $\blacksquare$

###  III-B Characterization of the follower’s best response

The following proposition establishes the characterization of the best
response $\mathrm{BR}_{2,\pi}$ in (3) for the follower and for a given
strategy $\gamma_{1}$ of the leader.

###### Proposition 2

Let $\gamma_{1}\in\Gamma_{1}$. The Best response of the follower is
characterized by:

$\begin{array}[]{l}(\mathcal{K}_{2,\pi}^{\mathrm{BR}}(\gamma_{1}),\widetilde{\mathcal{K}}_{2,\pi}^{\mathrm{BR}}(\gamma_{1}))\in\\\
\displaystyle\displaystyle\operatorname*{\text{argmax}\,}\limits_{\mathcal{K}_{2},\widetilde{\mathcal{K}}_{2}}\sum\limits_{k\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}p_{2k}\left(1-\frac{1}{\displaystyle\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)}\right)\\\
\displaystyle+\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\displaystyle\frac{\displaystyle\sqrt{\frac{\delta_{2k}}{p_{2k}}}\sum_{\ell\in\widetilde{\mathcal{K}}_{2}}\sqrt{p_{2\ell}\delta_{2\ell}}}{\displaystyle\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}\delta_{2\ell}\max\left(\displaystyle\frac{\gamma_{1\ell}}{\delta_{1\ell}}+\pi,1\right)}\right)\\\
\text{s.t. }\forall k\in\widetilde{\mathcal{K}}_{2},\\\
\frac{\displaystyle\sqrt{p_{2k}\delta_{2k}}\left[\displaystyle\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}\delta_{2\ell}\max\left(\displaystyle\frac{\gamma_{1\ell}}{\delta_{1\ell}}+\pi,1\right)\right]}{\displaystyle\sum_{\ell\in\widetilde{\mathcal{K}}_{2}}\sqrt{p_{2\ell}\delta_{2\ell}}}\\\
\hskip
113.81102pt>\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right).\end{array}$

Furthermore, any
$\gamma_{2,\pi}^{\mathrm{BR}}\in\mathrm{BR}_{2,\pi}(\gamma_{1})$ is defined
by: $\gamma_{2k,\pi}^{\mathrm{BR}}=$

$\left\\{\begin{array}[]{ll}0\hskip 102.43008pt\mbox{ if
}k\in\mathcal{K}\setminus{\mathcal{K}_{2,\pi}^{\mathrm{BR}}(\gamma_{1})},\\\
\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\hskip
14.22636pt\mbox{ if
}k\in\mathcal{K}_{2,\pi}^{\mathrm{BR}}(\gamma_{1})\setminus\widetilde{\mathcal{K}}_{2,\pi}^{\mathrm{BR}}(\gamma_{1}),\\\
\frac{\displaystyle\sqrt{p_{2k}\delta_{2k}}\left[\displaystyle\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2,\pi}^{\mathrm{BR}}(\gamma_{1})\setminus\widetilde{\mathcal{K}}_{2,\pi}^{\mathrm{BR}}(\gamma_{1})}\gamma_{2\ell}^{\mathrm{BR}}\right]}{\displaystyle\sum_{\ell\in\widetilde{\mathcal{K}}_{2,\pi}^{\mathrm{BR}}(\gamma_{1})}\sqrt{p_{2\ell}\delta_{2\ell}}},\
\mbox{ if }k\in\widetilde{\mathcal{K}}_{2,\pi}^{\mathrm{BR}}(\gamma_{1}).\\\
\end{array}\right.$ (6)

_$\mathbf{Proof.}$_ According to Proposition 1, the best response of Firm $2$
is characterized by the best part $\mathcal{K}_{2}\in 2^{\mathcal{K}}$ that
maximises

$\begin{array}[]{ll}&\displaystyle\max_{\gamma_{2}}\displaystyle\sum_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)\\\
&\text{s.t.
}\gamma_{2}\in\widetilde{\Gamma}_{2,\pi}(\mathcal{K}_{2},\gamma_{1}),\
\displaystyle\sum_{k=1}^{K}\gamma_{2k}\leq\mathrm{B_{2}}.\end{array}$ (7)

Let us exploit the KKT conditions by defining first the Lagrangian:
$\displaystyle\mathcal{L}(\gamma_{2},\mu_{2},\underline{\mu}_{2},\lambda_{2}):=\displaystyle\sum_{k\in\mathcal{K}\setminus\mathcal{K}_{2}}\underline{\mu}_{2k}\gamma_{2k}$

$\begin{array}[]{l}+\displaystyle\sum\limits_{k\in\mathcal{K}_{2}}p_{2k}\left(1-\frac{\delta_{2k}}{\gamma_{2k}}\right)-\lambda_{2}\left(\sum_{k=1}^{K}\gamma_{2k}-\mathrm{B_{2}}\right)\\\
\displaystyle+\sum_{k\in{\mathcal{K}}_{2}}\mu_{2k}\Bigg{(}\gamma_{2k}-\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\Bigg{)}\\\
\displaystyle-\sum_{k\in\mathcal{K}\setminus{\mathcal{K}_{2}}}\mu_{2k}\Bigg{(}\gamma_{2k}\displaystyle-\delta_{2k}\max\left(\displaystyle\displaystyle\frac{\gamma_{1k}}{\delta_{1k}},1\right)\Bigg{)}.\\\
\displaystyle\end{array}$

Let us denote by $\gamma_{2k,\pi}^{\star}$, $\mu_{2k,\pi}^{\star}$,
$\lambda_{2,\pi}^{\star}$ and $\underline{\mu}_{2k,\pi}^{\star}$ the variables
that verify the first-order optimality condition.

For all $k\in{\mathcal{K}}\setminus\mathcal{K}_{2}$,
$\displaystyle\frac{\partial\mathcal{L}}{\partial\gamma_{2k}}=-\lambda_{2,\pi}^{\star}-\mu_{2k,\pi}^{\star}+\underline{\mu}_{2k,\pi}^{\star}=0,$
then $\mu_{2k,\pi}^{\star}=0,\
\underline{\mu}_{2k,\pi}^{\star}=\lambda_{2,\pi}^{\star}>0\text{ and
}\gamma_{2k,\pi}^{\star}=0.$

For all $k\in\mathcal{K}_{2}$,
$\displaystyle\frac{\partial\mathcal{L}}{\partial\gamma_{2k}}=\displaystyle\frac{p_{2k}\delta_{2k}}{({\gamma}_{2k,\pi}^{\star})^{2}}-\lambda_{2,\pi}^{\star}-\mu_{2k,\pi}^{\star}=0$,
then
${\gamma}_{2k,\pi}^{\star}=\displaystyle\sqrt{\frac{p_{2k}\delta_{2k}}{(\lambda_{2,\pi}^{\star}+\mu_{2k,\pi}^{\star})}}.$
Let $\widetilde{\mathcal{K}}_{2}\in 2^{\mathcal{K}_{2}}$ (that may be empty)
such that, $\widetilde{\mathcal{K}}_{2}\in\\{\widetilde{\mathcal{K}}\in
2^{\mathcal{K}_{2}}:\ \forall k\in\widetilde{\mathcal{K}},\
\mu_{2k,\pi}^{\star}=0\\}.$ For all $k\in\widetilde{\mathcal{K}}_{2}$, the
first-order optimality condition is verified when,
$\displaystyle\gamma_{2k,\pi}^{\star}=\sqrt{\frac{p_{2k}\delta_{2k}}{\lambda_{2,\pi}^{\star}}}>\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)$
and
$\gamma_{2k,\pi}^{\star}=\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right),\
\forall k\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}.$ Since
$\lambda_{2,\pi}^{\star}>0$, it follows that,
$\sum_{k=1}^{K}\gamma_{2k,\pi}^{\star}=\mathrm{B_{2}}.$ Hence,
$\displaystyle\sum_{\ell\in\widetilde{\mathcal{K}}_{2}}\gamma_{2\ell,\pi}^{\star}=\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}\gamma_{2\ell,\pi}^{\star}$
$\Rightarrow$
$\sqrt{\lambda_{2,\pi}^{\star}}=\frac{\displaystyle\sum_{\ell\in\widetilde{\mathcal{K}}_{2}}\sqrt{p_{2\ell}\delta_{2\ell}}}{\displaystyle\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}\gamma_{2\ell}^{\star}}.$
Thus, the solution of (7) is characterized by $\gamma_{2k,\pi}^{\star}=$

$\left\\{\begin{array}[]{ll}0\hskip 113.81102pt\mbox{ if
}k\in\mathcal{K}\setminus{\mathcal{K}_{2}},\\\
\delta_{2k}\max\left(\displaystyle\frac{\gamma_{1k}}{\delta_{1k}}+\pi,1\right)\hskip
24.18501pt\mbox{ if
}k\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2},\\\
\frac{\displaystyle\sqrt{p_{2k}\delta_{2k}}\left[\displaystyle\mathrm{B_{2}}-\sum_{\ell\in\mathcal{K}_{2}\setminus\widetilde{\mathcal{K}}_{2}}\gamma_{2\ell,\pi}^{\star}\right]}{\displaystyle\sum_{\ell\in\widetilde{\mathcal{K}}_{2}}\sqrt{p_{2\ell}\delta_{2\ell}}}\
\mbox{ if }k\in\widetilde{\mathcal{K}}_{2}.\\\ \end{array}\right.$

Finally, the best response of the follower is obtained by a selection of
$\mathcal{K}_{2}$ and $\widetilde{\mathcal{K}}_{2}$ solving the discrete OP as
stated in Proposition 2. $\blacksquare$

###  III-C Proof of Theorem 1

In view of (3) and Proposition 1, the utility of the leader can be
reformulated such as

$\displaystyle u_{1,\pi}^{S}(\gamma_{1})=\sum_{\mathcal{K}_{1}\in
2^{\mathcal{K}}}\Big{[}\sum_{k\in\mathcal{K}_{1}}p_{1k}(1-\frac{\delta_{1k}}{\gamma_{1k}})\mathbb{1}_{\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})}(\gamma_{1})\Big{]}$,
and we recall that $\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})$ is defined
before Theorem 1. Thus, we identify the same utility structure as in the proof
of Proposition 1 and at the Stackelberg equilibria, it can be written as

$u_{1,\pi}^{S}(\gamma_{1,\pi}^{S})=\max_{\mathcal{K}_{1}\in
2^{\mathcal{K}}}\max_{\gamma_{1}\in\Gamma_{1}}\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})\vspace*{-1em}$
(8)

where $\widehat{S}_{\pi}(\mathcal{K}_{1},\gamma_{1})$ is defined in (4).
Concerning the existence of the Stackelberg strategy, the result is mainly
based on the existence of a feasible solution in the set of constraints for
the leader. Since $\mathcal{K}_{1}=\emptyset$ and $\gamma_{1}=(0,\ldots,0)$ is
in the set of constraints, we derive that the game $\mathcal{G}_{\pi}^{S}$ has
at least one Stackelberg equilibrium.

In order to compute numerically (8) with well known solvers for convex
optimization problems, let us prove that
$\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})$ is a convex set. Let
$\gamma_{1}^{a}\in\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})$ and
$\gamma_{1}^{b}\in\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1})$. Hence $\forall
k\in\mathcal{K}_{1}$ and $i\in\\{a,b\\}$,

$\hskip
48.36958pt\displaystyle\frac{\gamma_{1k}^{i}}{\delta_{1k}}\geq\max\left(\displaystyle\frac{\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1}^{i})}{\delta_{2k}},1\right).$

From the best response characterization of the follower in Proposition 2 one
has that $\forall k\in\mathcal{K}_{1}$,
$\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1}^{i})=0$. Hence $\forall
k\in\mathcal{K}_{1}$, $\gamma_{1k}^{i}\geq\delta_{1k}$. Let $\tau\in(0,1)$ and
$\gamma_{1k}^{b}\geq\tau\gamma_{1k}^{a}+(1-\tau)\gamma_{1k}^{b}\geq\gamma_{1k}^{a}.$
From the monotony of the best response of the follower w.r.t the action of the
leader it follows that, $\forall k\in\mathcal{K}_{1}$

$\hskip
48.36958pt\mathrm{BR}_{2k,\pi}^{-}(\tau\gamma_{1k}^{a}+(1-\tau)\gamma_{1k}^{b})=0.$
Hence,
$\tau\gamma_{1k}^{a}+(1-\tau)\gamma_{1k}^{b}\in\widehat{\Gamma}_{1,\pi}(\mathcal{K}_{1}).$
Finally, (8) is a strictly convex OP that can be solved with numerical solver
for convex OP. Finally, at the Stackelberg equilibrium this analysis
guarantees that $\forall k\in\mathcal{K}_{1,\pi}^{S}$ (region where the leader
invest), $\mathrm{BR}_{2k,\pi}^{-}(\gamma_{1,\pi}^{S})=0$, and in view of
Proposition 2 we derive that the follower’s strategy is given by
$\gamma_{2k,\pi}^{S}=\frac{\displaystyle\sqrt{p_{2k}\delta_{2k}}\displaystyle}{\displaystyle\sum_{\ell\in{\mathcal{K}}_{2,\pi}^{S}}\sqrt{p_{2\ell}\delta_{2\ell}}}\mathrm{B_{2}}$
if $k\in{\mathcal{K}}_{2,\pi}^{S}$ and $\gamma_{2k,\pi}^{S}=0$ otherwise,
where $\mathcal{K}_{2,\pi}^{S}$ is defined in Theorem 1.

##  IV Numerical performance analysis

In this section, we illustrate the solutions of the Stackelberg game when both
firms play pessimistically over a network of $K=5$ regions. The parameters of
the game are given by: $p_{1}=(p_{11},\ldots,p_{15})=(1,2,3,4,5)$;
$p_{2}=(p_{11},\ldots,p_{25})=(2,3,1,5,4)$;
$\delta_{1}=(\delta_{11},\ldots,\delta_{15})=10^{-1}\times(5,4,3,2,1)$;
$\delta_{2}=(\delta_{21},\ldots,\delta_{25})=10^{-1}\times(1,2,3,4,5)$ and the
entry price for the follower is fixed at $\pi=10^{-6}$.

Budgets | $(\gamma_{11,\pi}^{S},\gamma_{21,\pi}^{S})$ | $(\gamma_{12,\pi}^{S},\gamma_{22,\pi}^{S})$ | $(\gamma_{13,\pi}^{S},\gamma_{23,\pi}^{S})$
---|---|---|---
(0.6,0.6) | (0.2,0) | (0.4,0) | (0,0)
(0.6,0.6) | (0,0) | (0,0) | (0.6,0)
(0.6,5) | (0.6,0) | (0,1.298) | (0,1.377)
(5,0.6) | (0.5,0) | (1,0) | (1.5,0)
(5,5) | (0.833,0) | (1.666,0) | (2.5,0)
Budgets | $(\gamma_{14,\pi}^{S},\gamma_{24,\pi}^{S})$ | $(\gamma_{15,\pi}^{S},\gamma_{25,\pi}^{S})$
---|---|---
(0.6,0.6) | (0,0.335) | (0,0.264)
(0.6,0.6) | (0,0.335) | (0,0.264)
(0.6,5) | (0,1.298) | (0,1.025)
(5,0.6) | (2,0) | (0,0.6)
(5,5) | (0,2.79) | (0,2.20)
TABLE I:  Budget allocation at the Stackelberg Strategy for different value of
$(\mathrm{B}_{1},\mathrm{B}_{2})$

Fig. 1 represents the revenues of the two firms in the utility region, when
they both apply Stackelberg’s pessimistic strategy and for different budget
values. In the case of equal budgets $\mathrm{B_{1}}=\mathrm{B_{2}}$ the
follower has a higher revenue at the Stackelberg equilibrium due to the
difference in the churn rate in regions 4 and 5. When one budget is much
higher the corresponding Firm gets a larger revenue. We also observe that, at
the Stackelberg equilibrium, there is no case where
$\mathrm{B_{1}}\leq\mathrm{B}_{2}$ such that
$u_{1,\pi}^{S}(\gamma_{\pi}^{S})\geq u_{2,\pi}^{S}(\gamma_{\pi}^{S})$. The
TABLE. I shows the Stackelberg strategy allocation for the couples
$(\mathrm{B}_{1},\mathrm{B}_{2})$ highlighted in Fig. 1.

Figure 1:  Revenue of each firm (right) at the Stackelberg equilbrium for
values of $(\mathrm{B_{1}},\mathrm{B_{2}})$ shown on the left.

##  V Conclusion

We have formulated a Stackelberg duopoly game in which two firms compete for a
larger market share when services propagate according a viral model. We have
characterized analytically the corresponding Stackelberg strategy in the
pessimistic/optimistic setting. To obtain this result, we have described the
best response map of the follower and proved that the best response map of the
leader can be found by solving a convex OP.

To summarize, the paper provides a Stackelberg game to mathematically model a
decision-making problem in a economic competition framework. The solution
proposes strategic budget allocation across different regions in order to get
a larger market share. It is noteworthy that we consider both the case of
emerging companies (like startups) that need to cross a ”barrier to entry” on
the market and well established companies that want to preserve and enlarge
their market share. Basically the paper gives insights on the strategic
resource allocation with a priori given limitations in a competition setting.

## References

  * [1] S. Tu and S. Neumann.  A viral marketing-based model for opinion dynamics in online social networks.  In ACM Web Conference, pages 1570–1578, 2022.
  * [2] I.C. Morărescu, V.S Varma, L. Buşoniu, and S. Lasaulce.  Space-time budget allocation policy design for viral marketing.  Nonlinear Analysis: Hybrid Systems, 37:100899, 2020.
  * [3] J. Liu, P-E Paré, A. Nedić, C-Y Tang, C-L Beck, and T. Başar.  On a continuous-time multi-group bi-virus model with human awareness.  In 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pages 4124–4129. IEEE, 2017.
  * [4] J. Liu, P-E Paré, A. Nedić, C-Y Tang, C-L Beck, and T. Başar.  Analysis and control of a continuous-time bi-virus model.  IEEE Transactions on Automatic Control, 64(12):4891–4906, 2019\.
  * [5] R. Srinivasan.  Winner-Takes-All Dynamics, pages 199–215.  Springer Singapore, Singapore, 2021.
  * [6] B.A Prakash, A. Beutel, R. Rosenfeld, and C. Faloutsos.  Winner takes all: competing viruses or ideas on fair-play networks.  In Proceedings of the 21st international conference on World Wide Web, pages 1037–1046, 2012.
  * [7] V.S Varma, I.C. Morărescu, S. Lasaulce, and S. Martin.  Opinion dynamics aware marketing strategies in duopolies.  In 2017 IEEE 56th annual conference on decision and control (CDC), pages 3859–3864. IEEE, 2017.
  * [8] S. Gracy, I.C. Morărescu, V.S Varma, and P.E Pare.  Analysis and on/off lockdown control for time-varying sis epidemics with a shared resource.  In 2022 European Control Conference (ECC), pages 1660–1665, 2022\.
  * [9] V. Taynitskiy, E. Gubar, D. Fedyanin, I. Petrov, and Q. Zhu.  Optimal control of joint multi-virus infection and information spreading.  IFAC-PapersOnLine, 53(2):6650–6655, 2020.
  * [10] O. Lindamulage De Silva, S. Lasaulce, and I.C. Morărescu.  On the efficiency of decentralized epidemic management and application to Covid-19.  IEEE Control Systems Letters, 6:884–889, 2022.
  * [11] K. Kandhway and J. Kuri.  How to run a campaign: Optimal control of SIS and SIR information epidemics.  Applied Mathematics and Computation, 231:79–92, 2014.
  * [12] C. Nowzari, V.M Preciado, and G.J Pappas.  Optimal resource allocation for control of networked epidemic models.  IEEE Transactions on Control of Network Systems, 4(2):159–169, 2015\.
  * [13] G.S Zaric and M.L Brandeau.  Dynamic resource allocation for epidemic control in multiple populations.  Mathematical Medicine and Biology, 19(4):235–255, 2002.
  * [14] V.M Preciado, M. Zargham, C. Enyioha, A. Jadbabaie, and G. Pappas.  Optimal vaccine allocation to control epidemic outbreaks in arbitrary networks.  In 52nd IEEE conference on decision and control, pages 7486–7491. IEEE, 2013.
  * [15] V.S Varma, Y. Hayel, and I.C Morărescu.  A non-cooperative resource utilization game between two competing malware.  IEEE Control Systems Letters, 7:67–72, 2022.
  * [16] O. Lindamulage de Silva, V.S Varma, I.C Morărescu, and S. Lasaulce.  Optimal influence budget allocation for viral marketing using a multiple virus SIS model.  https://hal.science/hal-04010064, 2023.
  * [17] H. Sofia Rodrigues and Manuel J. F.  Can information be spread as a virus? viral marketing as epidemiological model.  Mathematical methods in the applied sciences, 39:4780–4786, 2016\.
  * [18] K. Kandhway and J. Kuri.  How to run a campaign: Optimal control of SIS and SIR information epidemics.  Applied Mathematics and Computation, 231:79–92, 2014.
  * [19] H. Sofia Rodrigues and Manuel J? F.  Viral marketing as epidemiological model.  arXiv preprint arXiv:1507.06986, 2015.
  * [20] H. Demsetz.  Barriers to entry.  The American Economic Review, 72:47–57, 1982.
  * [21] K. Basu and N. Singh.  Entry-deterrence in stackelberg perfect equilibria.  International Economic Review, 31:61–71, 1990.
  * [22] V. Doshi, S. Mallick, et al.  Competing epidemics on graphs-global convergence and coexistence.  In IEEE INFOCOM 2021-IEEE Conference on Computer Communications, pages 1–10. IEEE, 2021.
  * [23] N. Niknami and J. Wu.  Competitive influence maximisation model with monetary incentive.  International Journal of Parallel, Emergent and Distributed Systems, 37(6):680–695, 2022.
  * [24] Z. Xu, A. Khanafer, and T. Başar.  Competition over epidemic networks: Nash and stackelberg games.  In 2015 American Control Conference (ACC), pages 2063–2068. IEEE, 2015.
  * [25] G. Leitmann.  On generalized stackelberg strategies.  Journal of optimization theory and applications, 26:637–643, 1978\.