# Surrogate model for Bayesian optimal experimental design in chromatography

Jose Rodrigo Rojo-Garcia  Computational Engineering, School of Engineering
Science, Lappeenranta-Lahti University of Technology, Finland  Heikki Haario
Computational Engineering, School of Engineering Science, Lappeenranta-Lahti
University of Technology, Finland  Tapio Helin  Computational Engineering,
School of Engineering Science, Lappeenranta-Lahti University of Technology,
Finland  Tuomo Sainio  Separation Science, School of Engineering Science,
Lappeenranta-Lahti University of Technology, Finland

###### Abstract

We applied Bayesian Optimal Experimental Design (OED) in the estimation of
parameters involved in the Equilibrium Dispersive Model for chromatography
with two components with the Langmuir adsorption isotherm. The coefficients
estimated were Henry’s coefficients, the total absorption capacity and the
number of theoretical plates, while the design variables were the injection
time and the initial concentration. The Bayesian OED algorithm is based on
nested Monte Carlo estimation, which becomes computationally challenging due
to the simulation time of the PDE involved in the dispersive model. This
complication was relaxed by introducing a surrogate model based on Piecewise
Sparse Linear Interpolation. Using the surrogate model instead the original
reduces significantly the simulation time and it approximates the solution of
the PDE with high degree of accuracy. The estimation of the parameters over
strategical design points provided by OED reduces the uncertainty in the
estimation of parameters. Additionally, the Bayesian OED methodology indicates
no improvements when increasing the number of measurements in temporal nodes
above a threshold value.
Keywords: Chromatography, parameter estimation, optimal experimental design,
Bayesian statistics.

##  1 Introduction

Chromatography is a method to separate chemical compounds in complex mixtures
by using a column packed with a selective adsorbent material. It is utilized
in vastly different scales; the amounts that can be processed range from
nanograms in analytical chemistry to hundreds notref, or thousands of tons per
year in process industry. The separation is based on differences in
propagation velocities of the compounds as they are eluted with a solvent that
is pumped through the bed. Molecular level interactions with the adsorbent
lead to a local distribution equilibrium (adsorption equilibrium) between the
mobile fluid phase and the stationary adsorbent phase. Adsorption equilibrium
is thus characteristic for the solution and the adsorbent but, at constant
temperature, it depends on the solute concentrations only. The retention of
the compounds in the column can be explained and predicted if the adsorption
equilibrium isotherm is known.

The inability to predict isotherms a priori from first principles, especially
for complex systems, underscores the necessity for their experimental
determination. Several methods have been developed for determining the
isotherms experimentally [55]. The accuracy and utility of the experimental
data and the models used to correlate them are paramount for model-based
design, a burgeoning approach in optimizing chromatographic processes. In this
context, the inverse method [5] often gains preference among researchers,
particularly when an existing mathematical model of the chromatographic system
is available and can be numerically solved through established computational
techniques.

Methodologically, a series of chromatographic experiments are undertaken to
empirically identify key parameters, both for the adsorption isotherms and for
the non-idealities in fluid flow, termed dispersion. Predictions of a dynamic
simulation model are matched with the experimental data that is presented as a
time-series of concentration values obtained at the chromatographic column’s
outlet. Yet, this approach is not without its challenges. Typically, a single
experiment can require an hour, and real-time, high-accuracy analysis may not
be possible with online detectors. Consequently, researchers may have to rely
on off-line analyses that are even more time-consuming. Moreover, the errors
produced by the numerical scheme used to model the chromatographic system can
be significant [34]. On the other hand, simplified isotherm models with less
parameters have been observed to yield more reliable (albeit physically less
meaningful) parameters with the inverse method [43].

Variability in model parameters introduces uncertainties that propagate
through to the resultant design, posing questions about its optimality.
Therefore, a rigorous exploration of how these uncertainties might be
mitigated through experiment design is both intriguing and timely. In this
article, we aim to delve into these intricacies, exploring both the
theoretical and practical aspects of model-based chromatographic design, with
an emphasis on the role of experimentally determined isotherms and the
implications of parameter uncertainty.

Mathematical description of chromatography process is typically based on
systems of hyperbolic or parabolic PDEs. In this study, we focus on the
Equilibrium Dispersive Model (EDM), which comprises a nonlinear parabolic
equation with an adsorption term governed by the Langmuir equations. This
model depends on parameters such as the number of theoretical plates, the
Langmuir constants, and Henry’s coefficient [20]. Two common approaches for
the numerical solution of this problem are the finite volume methods (FVM) and
the discontinuos Galerkin method (DG) [31, 32, 33]. One of the FVM methods
suggested in [32] is the Koren scheme [39], which will be employed in our
study.

Classical strategies for the model calibration by the inverse method are based
on optimization algorithms (see e.g. [6, 59]). Here, we adopt the Bayesian
approach to quantify uncertainty in our calibration process. The Bayesian
paradigm involves representing the unknown parameters as a posterior
probability distribution, derived from combining prior information with the
likelihood density through the Bayes’ theorem [19]. Posterior probability
density is often intractable and advanced numerical algorithms are needed to
characterize it such as the Markov Chain Monte Carlo (MCMC) methods [53]. The
Bayes approach has been used in other contexts of chromatography to estimate
parameters in mathematical models, both PDE-based or otherwise (see e.g. [12,
24, 26, 40, 62, 64]). Specifically, for the EDM, we focus on estimating
parameters such as the Langmuir coefficients, total adsorption capacity, and
the number of theoretical plates.

Mathematical modelling of chromatography process involves variables that can
be controlled when planning the experiments, the so-called design variables.
The poor specification of the design variables can significantly increase the
uncertainty in the posterior measure, emphasizing the role of experimental
design. For example, the controllable laboratory variables can include the
concentration of the injected substance and the injection time and the
sampling instances.

Optimal experimental design (OED) originally emerged in the context of
frequentist statistics and expanded to Bayesian approach (see e.g. [16, 19,
49]). OED has an increasing role in chromatography in order to maximize the
efficiency of the experiment, and to reduce waste and cost [27]. Majority of
previous literature considers experimental design through empirical response
models, obtained by fitting the model to measurement data, see e.g. [7, 52].
In this work, we leverage the underlying physico-chemical model of the system
to produce more precise analysis of the design problem. Consequentially, our
approach requires substantially larger computational effort. We note that, to
our knowledge, OED utilizing a physico-chemical model for chromatography has
been considered previously only in the frequentist context with D-OED criteria
in [50].

Bayesian OED aims to recover design that maximizes a given utility averaged
over the Bayesian joint distribution. Various different utilities have been
proposed, with the most prevalent being alphabetic criteria, notably the A-
and D-optimality. The A-optimality is based in the minimization of the
posterior variance [1, 2, 3]. The D-optimality criteria is connected to
maximizing information gain when updating prior into posterior distribution
[45, 10]. In this work, we focus on D-optimality, also known as expected
information gain, which for linear Gaussian problems reduces to minimizing the
determinant of the posterior covariance matrix.

The computational challenge in Bayesian OED emerges from the need to evaluate
and optimize an expected utility expressed as a high-dimensional integral.
Evaluating the integrand requires evaluation of the likelihood density and,
therefore, the physico-chemical model, which is expensive. Specialized
numerical algorithms, typically based on Monte Carlo methods, are needed due
to the nested structure of the expectation, see [54] and subsequent work e.g.
in [28, 29, 30].

Minimizing computational effort in Bayesian OED is a subject of active
research. One avenue of research seeks to utilize a Gaussian approximation,
such as the Laplace approximation, for the posterior distribution. However,
common error bounds for this approach are expressed in the asymptotic limit of
repetitive measurements or small noise. Additionally, the approximation
requires the Jacobian of the forward mapping, which is not always available
(see, e.g. [9, 15, 46, 47, 48, 25]). Another line of work aims to approximate
the forward mapping or the likelihood density using a surrogate that is fast
to evaluate but maintains high accuracy. It’s worth noting that the stability
of such surrogate modelling was recently demonstrated in [18]. Among the
various surrogate approaches, we highlight Polynomial Chaos Expansion (PCE)
over sparse grids studied in [63].

In this study, the primary challenge with PCE-based methods lies in their
sensitivity to discontinuities or regions of high gradient in the forward
mapping, which can lead to artificial oscillations when smooth basis elements
are utilized [42]. Instead, we employ Piecewise Sparse Linear Interpolation
(PSLI), which provides localized approximation of the forward mapping while
weakening the influence of domain dimension on the approximation rate [14,
35].

###  1.1 Our contribution

This work contributes to the estimation of parameters in the area of
chromatography and the design of experiments in several ways. We consider the
situation where rough estimates for the model parameters are available, e.g.,
after some preliminary experiments. However, the uncertainty of the model does
not allow a reliable optimization of the individual sampling instants, only
the concentration of the injected material and the injection time are used as
the design variables. We develop a surrogate modeling approach that reduces
the computational time so that the Bayesian optimization of experiments
becomes feasible. Especially,

  * •

The surrogate model enables the optimization and sensitivity analyses of the
design variables as well as the MCMC sampling of the model parameter
posteriors. Additionally, the accuracy of the surrogate model is verified to
be high enough for satisfactory results.

  * •

We establish the minimum number of uniformly distributed sampling instants
that enables the parameters to be estimated with small uncertainty. Also, we
show that after a sufficient number of samples, using a high enough
concentration injected, no relevant improvement in the parameter estimation
can be achieved with increased experiments.

###  1.2 Structure of the paper

The paper is organized as follow. Section 2 describes the governing equations
of the EDM in chromatography, and the FVM algorithm chosen. In the same
section we introduce the concept of Bayesian inverse problems, the
mathematical formulation of the parameters estimation in the chromatography
model, and the sampling techniques that are used for inversion. In section 3
we formalize the definition of D-OED. After that we explain how the surrogate
model based in PSLI, and mathematical properties. Additionally we explain the
numerical algorithms for evaluate the utility function and their order order
of convergence. The section 4 we define the numerical values of the model
parameters and design variables, with which we create synthetic data for
different numbers of samples in the sensor. Finally the simulation is done
with the true model and the surrogate model, the results are compared in both
cases. In section 5 we conclude the results presented in the analysis.

##  2 Mathematical preliminaries

###  2.1 Dispersive model for chromatography

In a chromatography process, dispersive models are particularly valuable for
describing the behavior of components within complex mixtures as they travel
through a chromatographic column. In them, the system is characterized by a
set of components $c_{i}$, where $i$ ranges from 1 to $N_{c}$, representing
the different species or compounds of interest. These components are governed
by a set of non-linear hyperbolic partial differential equations (PDEs). Each
component $c_{i}$ is described by the following governing equation

$\frac{\partial c_{i}}{\partial t}+F\frac{\partial q_{i}}{\partial
t}+u\frac{\partial c_{i}}{\partial
z}=D_{\text{app}}\frac{\partial^{2}c_{i}}{\partial
z^{2}},\quad(z,t)\in(0,L)\times(0,T),$ (2.1)

where $F$ is the volumetric flow rate, $q_{i}$ models the isotherms of the
system, $u$ is the linear velocity of the mobile phase and $D_{app}$ is the
axial dispersion coefficient. The equation (2.1) is subjected to the zero
initial conditions

$c_{i}(z,0)=0,\quad q_{i}(z,0)=0,\quad z\in(0,L).$ (2.2)

In particular, the latter initial condition indicates that at the column
inlet, no adsorption has occurred. A common model, connecting the isotherms
$q_{i}$ to the concentrations $c_{i}$, is given by the Langmuir equation

$q_{i}(z,t)=\frac{Q_{s}b_{i}c_{i}(z,t)}{1+\sum_{l=1}^{N_{c}}b_{l}c_{l}(z,t)},\quad(z,t)\in(0,L)\times(0,T),$
(2.3)

where $q_{i}$ represents the maximum adsorption capacity, $b_{i}$ is the
Langmuir constant for component $i$, and $Q_{s}$ is the total adsorption
capacity of the stationary phase.

Boundary conditions describing the inflow and outflow of fluids are specified
as follows

$\displaystyle c_{i}(0,t)$
$\displaystyle=c_{i,0}(t)+\frac{D_{i}}{u}\frac{\partial c_{i}}{\partial
z}(0,t),\quad i=1,2,\ldots,N_{c},$ $\displaystyle\frac{\partial
c_{i}}{\partial z}(L,t)$ $\displaystyle=0.$

Notice that at the column’s outlet ($z=L$), the spatial gradient of each
component is set to zero, ensuring no mass transfer at the exit. In what
follows, we consider injections $c_{i,0}(t)$ satisfying

$c_{i,0}(t)=\begin{cases}c_{i}^{\text{Feed}}&\text{for }t\leq
t^{\text{inj}},\\\ 0&\text{for }t>t^{\text{inj}},\end{cases}\quad
i=1,\ldots,N_{c}.$ (2.4)

A common assumption is the dependency between the dispersion parameter
$D_{app}$ and the velocity $u$, and a common choice is the linear relation
[20]

$D_{\mathrm{app}}=\frac{Lu}{2N_{tp}},$

here the parameter $N_{tp}$ refers to theoretical plates formed during the
process.

For numerical purposes, we rescale spatial and time variables according to

$t=\frac{L}{u}\tau,\quad\tau\in[0,\Upsilon)\quad\text{and}\quad z=Ly,\quad
y\in[0,1],$

where $\Upsilon=\frac{uT}{L}$. Then, the problem is transformed to

$\frac{\partial c_{i}}{\partial\tau}+F\frac{\partial
q_{i}}{\partial\tau}+\frac{\partial c_{i}}{\partial
y}=\frac{1}{2N_{tp}}\frac{\partial^{2}c_{i}}{\partial y^{2}},$ (2.5)

with initial conditions

$c_{i}(0,y)=0,\quad q_{i}(0,y)=0,$ (2.6)

and boundary conditions

$\left.c_{i}\right|_{y=0}=c_{i,0}+\left.\frac{1}{2N_{tp}}\frac{\partial
c_{i}}{\partial y}\right|_{y=0},\quad\frac{\partial c_{i}(1,\tau)}{\partial
y}=0,$ (2.7)

and where

$c_{i,0}(\tau)=\left\\{\begin{array}[]{ll}c_{i}^{\text{Feed }}&\text{ for
}\tau\leq\tau^{\mathrm{inj}}\\\ 0&\text{ for
}\tau>\tau^{\mathrm{inj}}\end{array}\right.$ (2.8)

With this change of variables, the transformed model depends on $u$ only
through the new time variable $\tau$ and it simplifies the structure of the
model.

######  Remark 1.

The regularity of the solutions the PDE system depends on the regularity of
the boundary conditions and the initial condition: for regular enough
conditions the solution is differentiable. More mathematically (see [41, 44,
65]), fixing the model parameters and assuming
$c_{i}(0,\cdot)\in(L^{2}(0,L))^{N_{c}}$ and $c_{i,0}\in(L^{2}(0,T))^{N_{c}}$,
there exist a unique and stable solution in
$c_{i}\in(W_{2}^{2,1}((0,L)\times(0,T)))^{N_{c}}$.

###  2.2 Numerical implementation by the Koren scheme

The Koren scheme is a FVM method based on a piecewise-polynomial interpolation
for the flux, and a flux limiter proposed by Sweby [58]. Consider regular grid
$y_{m}$, $m=0,...,M$ on a unit interval with $y_{0}=0$ and $y_{M}=1$. We
denote the mid-points of the grid by $y_{m+1/2}$ for $m=0,...,M-1$, and the
stepsize by $\Delta y$, respectively. For what follows, let us also denote
$\bar{c}=(c_{i})_{i=1}^{N_{c}}:(0,1)\times(0,\Upsilon)\to\mathbb{R}^{N_{c}}$,
and the flux as $f(\bar{c})=\bar{c}.$ By chain rule, we have

$\frac{\partial
q_{i}}{\partial\tau}(y,\tau)=F\sum_{j=1}^{N_{c}}\frac{dq_{i}}{dc_{j}}(y,\tau)\frac{\partial
c_{j}}{\partial\tau}(y,\tau)=F(\nabla
q_{i})(y,\tau)\cdot\frac{\partial\bar{c}}{\partial\tau}(y,\tau)$

with the convention $(\nabla
q_{i})=(\frac{dq_{i}}{dc_{j}})_{j=1}^{N_{c}}:(0,1)\times(0,\Upsilon)\to\mathbb{R}^{N_{c}}$.
In consequence, the vectorized form of the system (2.5) is

$(I+FQ)\frac{\partial\bar{c}}{\partial\tau}=-\frac{\partial\bar{c}}{\partial
y}+\frac{1}{2N_{tp}}\frac{\partial^{2}\bar{c}}{\partial y^{2}},$

where $I$ stands for an $N_{c}\times N_{c}$ identity matrix and $Q$ is
composed of the elements

$Q_{ij}=\frac{dq_{i}}{dc_{j}}:(0,1)\times(0,\Upsilon)\to\mathbb{R},\quad 1\leq
i,j\leq N_{c}.$

The Koren scheme is then given by

$\frac{d\bar{c}}{d\tau}(y_{m},\tau)\\\
=-\left(I+FQ(y_{m},\tau)\right)^{-1}\left[\frac{f_{m+\frac{1}{2}}(\tau)-f_{m-\frac{1}{2}}(\tau)}{\Delta
y}-\frac{1}{2N_{tp}\Delta y}\left(\frac{\partial\bar{c}}{\partial
y}(y_{m+\frac{1}{2}},\tau)-\frac{\partial\bar{c}}{\partial
y}(y_{m-\frac{1}{2}},\tau)\right)\right],$ (2.9)

where

$f_{\frac{1}{2}}(\tau)=\bar{c}(0,\tau),\quad
f_{\frac{3}{2}}(\tau)=\bar{c}(y_{1},\tau),\qquad
f_{N+\frac{1}{2}}(\tau)=\bar{c}(y_{N},\tau),$

and

$f_{m+\frac{1}{2}}=f_{m}+\frac{1}{2}\phi\left(r_{m+\frac{1}{2}}\right)\left(f_{m}-f_{m-1}\right).$

Above, the term $\phi$ is the flux limiter and it is defined as

$\phi\left(r_{m+\frac{1}{2}}\right)=\max\left(0,\min\left(2r_{m+\frac{1}{2}},\min\left(\frac{1}{3}+\frac{2}{3}r_{m+\frac{1}{2}},2\right)\right)\right),$

where

$r_{m+\frac{1}{2}}=\frac{f_{m+1}-f_{m}+\eta}{f_{m}-f_{m-1}+\eta},$

represents the ratio between consecutive fluxes with some small $\eta>0$ avoid
numerical singularities. Furthermore, we utilize the central approximation

$\frac{\partial\bar{c}}{\partial
y}(y_{m+\frac{1}{2}},\tau)\approx\frac{\bar{c}(y_{m+1},\tau)-\bar{c}(y_{m},\tau)}{\Delta
y}.$

Consider the local truncation error defined by

$\bar{E}(\tau):=\frac{d\bar{c}}{d\tau}(y_{m},\tau)\\\
+\left(I+FQ(y_{m},\tau)\right)^{-1}\left[\frac{f_{m+\frac{1}{2}}(\tau)-f_{m-\frac{1}{2}}(\tau)}{\Delta
y}-\frac{1}{2N_{tp}\Delta y}\left(\left(\frac{\partial\bar{c}}{\partial
y}\right)_{m+\frac{1}{2}}-\left(\frac{\partial\bar{c}}{\partial
y}\right)_{m-\frac{1}{2}}\right)\right].$

The Koren scheme has spatial consistency order 2, i.e.,
$\|\bar{E}(\tau)\|\leq\mathcal{O}(|\Delta y|^{2})$ (see e.g. [31, 33, 39]).
Following the methodology in [31], we use Runge-Kutta 4 (RK4) for temporal
discretization due it is well known that RK4 has consistency order 5.

###  2.3 Bayesian inference

Let us now formalize our experiment. The solution $c_{i}$, $i=1,...,N_{c}$ of
the system (2.1)-(2.8) is specified by the values of parameters $F$, $N_{tp}$,
$b_{1},...,b_{N_{c}}$, $Q_{s}$, $t^{\mathrm{inj}}$,
$c_{1}^{\mathrm{Feed}},...,c_{N_{c}}^{\mathrm{Feed}}$. For what follows, let
${\mathcal{P}}$ stand for the Cartesian product of all individual parameter
domains. We next decompose ${\mathcal{P}}$ into three components
${\mathcal{P}}=\Theta\times\mathcal{D}\times\mathcal{E}$, where $\Theta$,
$\mathcal{D}$ and $\mathcal{E}$ denote the sets of parameters of interest,
design parameters in our experiment and fixed variables, respectively. The
isotherm parameters $b_{j}$, $Q_{s}$ and $N_{tp}$ are the components of
$\Theta$ that need to be estimated. The design variables $\mathcal{D}$ are the
feed concentrations $c_{j}^{\mathrm{Feed}}$ and the injection time
$t^{\mathrm{inj}}$. And finally, the only fixed variable in the model is $F$.
Below, the process of estimating these parameters of interest is referred to
as the inverse problem. We define a mapping from the model parameters and
design variables to the model solutions as
$G:\Theta\times\mathcal{D}\rightarrow(W_{2}^{2,1}((0,1)\times(0,\Upsilon)))^{N_{c}}$
as the unique solution of the system (2.5)-(2.8) specified by Remark 1. In
addition, let
${\mathcal{O}}(\cdot):(W_{2}^{2,1}((0,1)\times(0,\Upsilon)))^{N_{c}}\rightarrow\mathbb{R}^{K}$
denote the observation operator evaluating the solution at a grid of $K$ time
points on $(0,1)\times(0,\Upsilon)$. We remark that in chromatography we only
do measurements at the outlet $y=1$.

In the experiment, we wish to determine the unknown parameters
$\theta\in\Theta$ from the noisy observations

${\bf c}=\mathcal{O}(G(\theta;d))+\eta=\mathcal{G}(\theta;d)+\eta,$ (2.10)

where ${\bf c},\eta\in\mathbb{R}^{K}$ represent vectors of point evaluations
of the concentrations and noise, respectively, and we write
$\mathcal{G}=\mathcal{O}\circ G:\Theta\times\mathcal{D}\to\mathbb{R}^{K}$. The
realization of the noise process is assumed to be drawn from the Gaussian
distribution ${\mathcal{N}}(0,\sigma^{2}I)$, i.e., each point evaluation is
independent with the same noise level $\sigma$. The value of $\sigma$ is
supposed to be known, as estimated by fitted residuals.

Due to the presence of noise in the data, there is inherent uncertainty
associated with any estimator of $\theta$. The Bayesian inference provides a
paradigm to assess this uncertainty by updating any prior distribution of
$\theta$ into a posterior distribution. We assume here Gaussian measurement
noise. The parameter uncertainty is then formalized by the next theorem.

######  Theorem 2 (Bayes’ theorem).

Let $p_{0}(\theta)$ denote the prior probability density of $\theta\in\Theta$.
Then the posterior density function $p(\theta|{\bf c};d)$ (the density of
$\theta$ given ${\bf c}\in\mathbb{R}^{K}$) corresponding to the measurement
(2.10) is given by

$p(\theta\mid{\bf c},d)=\frac{1}{Z({\bf
c})}\exp\left(-\frac{\left\|\mathcal{G}(\theta;d)-{\bf
c}\right\|^{2}}{2\sigma^{2}}\right)p_{0}(\theta),$

where

$Z({\bf c})=\int_{\Theta}\exp\left(-\frac{\left\|\mathcal{G}(\theta;d)-{\bf
c}\right\|^{2}}{2\sigma^{2}}\right)p_{0}(\theta)d\theta.$

To explore the posterior posterior distribution, Markov Chain Monte Carlo
Methods (MCMC) are commonly employed, as they can generate samples from the
posterior distribution without the need for the normalization constant, as
discussed in [53]. Producing a sample size that accurately represents the
final distribution may demand substantial computational effort. Various
adaptive techniques have been developed, as discussed in [17, 22, 23]. One
pragmatic and efficient approach is the DRAM algorithm, as detailed in [21],
which is supported by a MATLAB toolbox and has demonstrated satisfactory
results.

##  3 Optimal experimental design

In this section, we formulate Bayesian approach to experimental design. The
uncertainty variables are $\theta=(b_{1},b_{2},Q_{s},N_{tp})$, the design
variables are $d=(\tau^{inj},c^{Feed})$ and the measurements are given by a
vector of size $K=2N_{s}$ with evaluations of $c_{1}$ and $c_{2}$ in the
outlet $y=1$ and at times $\tau_{1},\tau_{2},\cdots,\tau_{N_{s}}$ as

${\bf
c}=(c_{1}(1,\tau_{1}),c_{1}(1,\tau_{2}),\cdots,c_{1}(1,\tau_{N_{s}}),c_{2}(1,\tau_{1}),c_{2}(1,\tau_{2}),\cdots,c_{2}(1,\tau_{N_{s}})).$

Bayesian optimal experimental design aims to maximize the expected utility
associated with a given design variable [16]. In more precise mathematical
notation, the objective is to maximize the utility function

$U(d)=\int_{\mathbb{R}^{K}\times\Theta}u(\theta,{\bf c};d)p(\theta,{\bf
c};d),$

where $p(\theta,{\bf c};d)$ is the Bayesian joint density function of $\theta$
and ${\bf c}$ on $\mathbb{R}^{K}\times\Theta$ and $u$ is the information gain.

Here, we consider the D-OED criteria based in the Kullback–Leibler divergence
(KL) between the posterior and prior measures, i.e

$u(\theta,{\bf c};d)=D_{KL}(p(\cdot\mid{\bf c};d)\;\|\;p_{0}),$

and the expected utility function is then reduced as

$U(d)=\mathbb{E}D_{KL}(p(\cdot\mid{\bf
c};d)\;\|\;p_{0})=\iint_{\mathbb{R}^{K}\times\Theta}p(\theta|{\bf
c};d)\log\left(\dfrac{p(\theta|{\bf c};d)}{p_{0}(\theta)}\right)d\theta p({\bf
c};d)d{\bf c},$ (3.1)

where $p({\bf c};d)$ is the marginal density of the measurement ${\bf c}$.
This quantity indicates how similar are posterior and prior measures in the
design node $d$, and a big value implies a more informative posterior
distribution in that node.

A main challenge of Bayesian experimental design is the numerical
approximation of the above double integral. For linear inverse problems with
Gaussian prior distribution, the expected information gain (3.1) can be
evaluated in closed form and it is reduced to the computation of the log-
determinant of the posterior covariance matrix. In the inverse problem
literature, this property has been widely utilized in various contexts; see,
e.g., [1, 4, 57].

For non-linear forward mappings that do not have sufficient linearization
properties, such as our problem here, the integral in (3.1) needs to be
evaluated approximatively. While Monte Carlo based algorithms are a natural
choice, they tend to require high computational effort due to the inherent
nested sampling task, which we discuss more below. Before considering
evaluation, let us note that the computational cost is dominated by the effort
needed to evaluate the forward mapping. Well-design surrogate models can
reduce this cost.

###  3.1 Double-loop Monte Carlo integration

The expected utility in (3.1) can be rephrased as

$U(d)=\int_{\mathbb{R}^{K}\times\Theta}\\{\log[p({\bf
c}|\theta;d)]-\log[p({\bf c};d)]\\}p({\bf c},\theta;d)d\theta d{\bf c},$ (3.2)

where $p({\bf c};d)$ refers to the marginal density of the observation, i.e.,
the evidence. The computational challenge in estimating (3.2) is that the
integrand function intractable (excluding special cases) and cannot be
directly approximated by conventional Monte Carlo methods. In literature, the
traditional approach is to apply double-loop (i.e. nested) Monte Carlo
integration as follows: first, an ensemble $\\{(\theta^{k},{\bf
c}^{k})\\}_{k=1}^{M}$ is generated from the joint distribution. Second, noting
the identity

$p\left({\bf c}^{k};d\right)=\int_{\Theta}p\left({\bf
c}^{k}|\theta;d\right)p_{0}(\theta)d\theta,$

the evidence density $p({\bf c}^{k};d)$ is approximated using a prior ensemble
$\\{\theta^{j,k}\\}_{j=1}^{J}$ for each $k=1,...,M$. This leads to the nested
MC approximation

$U_{M}(d)=\frac{1}{M}\sum_{k=1}^{M}\log\left[p\left({\bf
c}^{k}|\theta^{k};d\right)\right]-\frac{1}{JM}\sum_{k=1}^{M}\sum_{j=1}^{J}\log\left[p\left({\bf
c}^{k}|\theta^{j,k};d\right)\right].$

Notice that the estimator $U_{M}$ is biased positively, i.e.
$\mathbb{E}[U_{M}(d)-U(d)]\geq 0$ and the asymptotic mean squared error (MSE)
is of the order $\mathbb{E}|U_{M}(d)-U(d)|^{2}=\mathcal{O}(1/M+1/J^{2})$ (see
[54]). By scaling $J\propto\sqrt{M}$, we obtain MSE error ${\mathcal{O}}(1/M)$
with the computational cost of $M^{3/2}$ number of likelihood evaluations and,
therefore, evaluations of the forward operator ${\mathcal{G}}$. The line of
work [28, 29, 30] suggest replacing the nested ensemble
$\\{\theta^{j,k}\\}_{j=1}^{J}$ for each $k=1,...,M$, by the original prior
ensemble $\\{\theta^{k}\\}_{k=1}^{M}$, reducing the computational effort. The
resulting method is summarized in Algorithm 1.

Algorithm 1 Accelerated double-loop Monte Carlo estimation of expected utility
following [28]

Ensemble $\\{\theta_{k}\\}_{k=1}^{M}\sim p_{0}(\theta)$ i.i.d. Design variable
$d$.

$U(d)$

$U(d)\leftarrow 0$

for $k=1:M$ do

$\widehat{{\bf c}}^{k}\leftarrow{\mathcal{G}}(\theta^{k};d)$ $\triangleright$
Map all elements of the ensemble.

end for

for $k=1:M$ do $\triangleright$ Outer loop

${\bf c}^{k}\sim\mathcal{N}(\widehat{{\bf c}}^{k},\sigma^{2}I)$

$\log\left[p\left({\bf c}^{k}\mid\
\theta^{k},d\right)\right]\leftarrow-\dfrac{\|{\bf c}^{k}-\widehat{{\bf
c}}^{k}\|_{2}^{2}}{2\sigma^{2}}$

$p\left({\bf c}^{k}\mid d\right)\leftarrow 0$

for $j=1:M$ do $\triangleright$ Inner loop

$p\left({\bf c}^{k}\mid d\right)\leftarrow p\left({\bf c}^{k}\mid
d\right)-\frac{1}{M}\exp\left(-\dfrac{\|{\bf c}^{k}-\widehat{{\bf
c}}^{j}\|_{2}^{2}}{2\sigma^{2}}\right)$

end for

$U(d)\leftarrow U(d)+\frac{1}{M}\left(\log\left[p\left({\bf c}^{k}\mid\
\theta^{k},d\right)\right]-\log\left[p\left({\bf c}^{k}\mid
d\right)\right]\right)$

end for

###  3.2 Surrogated model based in Piecewise Linear Interpolation

In our context, the main computational effort associated to the double-loop MC
algorithm 1 arises from the repeated evaluations of the likelihood density
and, therefore, of the forward operator ${\mathcal{G}}$. Consequently, it is
tempting to replace ${\mathcal{G}}$ with a surrogate that offers both accurate
approximation and reduced computational cost. For computational purposes, we
consider the approximative observational model (2.10) by an approximation

${\bf c}_{N}={\mathcal{G}}_{N}(\theta;d)+\eta$

with a surrogate mapping
${\mathcal{G}}_{N}:\Theta\times\mathcal{D}\rightarrow\mathbb{R}^{K}$. The
parameter $N$ identifies the approximation rate, which here will correspond to
the discretization grid. In [18], it was demonstrated that the EIG (3.1) is
stable w.r.t. perturbations in the likelihood. In particular, it was shown
that the induced error in EIG evaluation is bounded by the squared
approximation error ${\mathcal{G}}(\theta;d)-{\mathcal{G}}_{N}(\theta;d)$
averaged over the prior.

In this work, we utilize piecewise sparse linear interpolation (PSLI) [14],
which is well-suited for surrogate modelling of mappings with limited
smoothness. For their rigorous construction, let us begin by introducing a
one-dimensional piecewise linear interpolant $U[f]$ of a function
$f:[-1,1]\rightarrow\mathbb{R}$. This interpolant is defined over a set of
nodes $-1=x_{0}<x_{1}<x_{2}<\cdots<x_{m}=1$ as follows:

$U[f](x)=\sum_{j=0}^{m}f(x_{j})a_{j}(x),\quad x\in[-1,1],$

where the basis is given by the hat functions

$a_{j}(x)=\begin{cases}1-\frac{n\left|x-x_{j}\right|}{2},&\text{ if
}\left|x-x_{j}\right|<\frac{2}{n}\\\ 0,&\text{ otherwise. }\end{cases}$

In what follows, we utilize the Clenshaw–Curtis quadrature rule defined by
nodes $x_{j}=-\cos(j\pi/m)$. Generalizing to higher dimensions, we denote
$U^{i}$ for the one-dimensional interpolant acting on the $i^{th}$ coordinate.
The multivariate interpolant on $[-1,1]^{r}$ is then constructed as a tensor
product

$(U^{i_{1}}\otimes...\otimes
U^{i_{r}})[f](x^{1},\cdots,x^{r})=\sum_{j_{1}=0}^{m_{1}}\sum_{j_{2}=0}^{m_{2}}\cdots\sum_{j_{r}=0}^{m_{r}}f(x_{j_{1}}^{1},\cdots,x_{j_{r}}^{r})(a_{j_{1}}(x^{1})\cdots
a_{j_{r}}(x^{r})),$ (3.3)

where $(x^{1},x^{2},\cdots,x^{r})\in[-1,1]^{r}$.

The drawback of (3.3) is the exponentially growing number of function
evaluations w.r.t. the dimension $r$ when $m_{1}=...=m_{r}$. Sparse grid
methods [14] provide a remedy for this curse of dimension by limiting the sum
in (3.3) to a carefully designed sparse subset of the grid points without
sacrificing too much accuracy. The Smolyak formulas [56] construct the sparse
grid leading to a hierarchical basis $\\{a_{j_{k}}\\}_{k=1}^{r}$ and
preserving the interpolation properties of $r=1$ when extended to higher
dimensions.

With the convention $U^{0}=0$, the Smolyak formula is defined as follows:
given $q\geq r$ we set

$\displaystyle S_{q,r}$ $\displaystyle=$ $\displaystyle\sum_{|\mathbf{i}|\leq
q}\left((U^{i_{1}}-U^{i_{1}-1})\otimes\cdots\otimes(U^{i_{r}}-U^{i_{r}-1})\right)$
(3.6) $\displaystyle=$ $\displaystyle\sum_{q-r+1\leq|\mathbf{i}|\leq
r}(-1)^{r-|\mathbf{i}|}\cdot\left(\begin{array}[]{c}r-1\\\
q-|\mathbf{i}|\end{array}\right)\cdot\left(U^{i_{1}}\otimes\cdots\otimes
U^{i_{r}}\right),$

where ${\bf i}=(i_{1},...,i_{r})\in{\mathbb{N}}^{r}$ and $|{\bf
i}|=i_{1}+...+i_{r}$. For the second identity, see e.g. [60, Lemma 1]. The
resulting grid consists of $N(q,r)\leq 2^{q}\frac{(q-1)^{r-1}}{(r-1)!}$ points
(see e.g. [51]); the arrangement is illustrated on a two-dimensional domain in
Figure 1.

Figure 1: Smolyak-Sparse Grids in 2D. The number of nodes are $N=65,145,321$
for left, center and right respectively.

The interpolant $S_{q,r}$ has well-established approximation properties. For
the characretization, we introduce the function space (see [13] and [51])

$F_{d}^{k}=\begin{cases}C^{k}([-1,1])&\text{ if }r=1\\\
\left\\{f:[-1,1]^{r}\rightarrow\mathbb{R}\mid D^{\alpha}f\text{ continuous if
}\alpha_{i}\leq k\text{ for all }i\right\\}&\text{ if }r>1\end{cases}$

######  Theorem 3 ([8, Remark 10]).

Given $f\in F_{r}^{k}$, $k\in\\{1,2\\}$, the interpolant $S_{q,r}$ satisfies

$\left\|f-S_{q,r}f\right\|_{\infty}\leq c_{r,k}N^{-k}(\log
N)^{(k+1)(r-1)}\left\|f\right\|_{C^{k}},$

where $N=N(q,r)$ is the number of grid points evaluated in $S_{q,r}$.

We note that piecewise linear interpolation is not order optimal for $k>2$
and, therefore, the rate in Theorem 3 does not extend beyond $k=2$, see [8,
Remark 10]. An efficient algorithmic implementation is discussed in [35, 36],
and a free toolbox is available [37, 38]. Moreover, Theorem 3 does not
guarantee convergence for continuous functions. Nonetheless, continuity is
typically the minimum requirement for interpolation, and numerical experiments
in [35, 51] show convergence for continuous benchmark functions. Finally, we
define our surrogate model by constructing

${\mathcal{G}}_{N}(\cdot;d_{j})=(S_{q,r}{\mathcal{G}}_{k}(\cdot;d_{j}))_{k=1}^{K},$

with ${\mathcal{G}}=({\mathcal{G}}_{k})_{k=1}^{K}$, on a grid of points
$d_{j}\in{\mathcal{D}}$, where $N=N(q,r)$. The full mapping for all
$d\in{\mathcal{D}}$ is obtained by linear interpolation. We note that while
theoretically establishing smoothness of the forward mapping ${\mathcal{G}}$
is outside the scope of this work, no issues were encountered in our numerical
simulations.

##  4 Numerical simulations

###  4.1 Problem formulation

In our numerical simulations, we model a chromatography process for two
components, i.e., $N_{c}=2$. Moreover, we assume that our measurements are
obtained on an equidistant temporal grid with $N_{s}$ nodes and the
measurement data is contaminated with a normally distributed noise with
standard deviation $\sigma$.

We utilize a conservative prior model composed of independent uniform
distributions for every component of the parameter $\theta$. That is, for
every parameter $\theta_{i}$, we define a uniform prior in each interval
$\left[\alpha_{i},\beta_{i}\right]\subset\mathbb{R}_{+}$, i.e.,
$\theta_{i}\sim\mathcal{U}\left(\left[\alpha_{i},\beta_{i}\right]\right)$.
Furthermore, we can restrict the domain to the Cartesian product of the
intervals $\Theta=\bigtimes_{i=1}^{4}\left[\alpha_{i},\beta_{i}\right]$. We
assume that these intervals can be reliably specified from independent
information, e.g. a rough preliminary experiment. We note that similar
methodology was utilized in [61].

As earlier, our design variables consist of the injection time and
concentration. Their range is limited due to the experimental conditions at
the laboratory. In consequence, we specify the design domain as the Cartesian
product
$\mathcal{D}=\left[\tau_{0}^{\text{inj}},\tau_{F}^{\text{inj}}\right]\times\left[c_{0}^{\text{Feed}},c_{F}^{\text{Feed}}\right]$.

###  4.2 Simulations

For this study, we exclusively utilized synthetic data, which we divided into
three experiments. In each experiment, measurement data was generated on an
equidistant grid with $N_{s}=8,15$, and $20$ temporal nodes within the time
interval $[0.5,9.5]$.

In the spirit of [11], we specified the true data-generating parameters as
$b_{1}=0.05$, $b_{2}=0.10$, $Q_{s}=10$. Moreover, we set $N_{tp}=70$. The
standard deviation of noise was $\sigma=0.05$, and other parameters satisfied
$\Upsilon=10$ and the constant $F=1.5$. The intervals of the prior
distribution and for the design space are shown in the tables 1 and 2,
respectively.

$Parameter$ | Lower Bound | Upper Bound | Real parameter
---|---|---|---
$b_{1}$ | 0.02 | 0.08 | 0.05
$b_{2}$ | 0.03 | 0.17 | 0.1
$Q_{s}$ | 8 | 11 | 10
$N_{tp}$ | 50 | 180 | 70
Table 1: Uncertainty Parameters. $Parameter$ | Lower Bound | Upper Bound
---|---|---
$\tau^{inj}$ | 0.05 | 3
$c^{Feed}$ | 1 | 15
Table 2: Design Parameters.

In all cases we estimated the utility function with a high-resolution forward
mapping and with PSLI for $N=1105$ training nodes. Every training process was
carried out on an equidistant grid over the design space with 14 nodes in each
direction. The average times needed for train the surrogate model is given in
Table 3. We compare the times invested in training the surrogate model by
using the high resolution model. In the case of the utility function, we can
see how the training time of the surrogate is one order of magnitude smaller
than the brute-force evaluation with real model. The time needed to calculate
a parameter chain of length 80 000 with MCMC is three orders of magnitude
smaller than that using the real model. The utility functions are presented in
the Figure 5.

In Figure 5, we observe periodic behavior for $N_{s}=8$, which gradually
diminishes for higher number of temporal nodes. This suggests that the
likelihood effectively captures the dominant components of the model relative
to the sensor samples. As expected, fewer sensor samples require precise
positioning to optimize information capture. Furthermore, there is a
noticeable increase in the utility function with respect to the concentration
$c^{\text{Feed}}$ across all cases. Additionally, the behavior concerning
$\tau^{\text{inj}}$ appears to become somewhat independent for $N_{s}=15,20$.

To validate the utility function results, we selected 6 design nodes (L, M, N,
R, T, and X) and computed MCMC samples using both the true model and the
surrogate model at these points. These nodes were chosen to facilitate
comparisons across rows and columns, demonstrating that the optimal posterior
distribution does not always directly correlate with $\tau^{inj}$ and
$c^{Feed}$. The posterior sampling was carried out with the DRAM algorithm
[21], comprising 80 000 simulations with the first 30 000 samples discarded to
mitigate the burn-in effect. The results are illustrated in figures 6, 7, and
8.

For $N_{s}=8$, we observe periodic behavior in the utility function concerning
the time of injection. The posterior distribution with less uncertainty tends
to concentrate around point $X$, whereas the distribution with more
uncertainty is typically centered at design node $M$. Additionally, we
anticipate lower uncertainty in the posterior distribution at point $N$
compared to $L$, and at $R$ compared to $T$.

The periodicity is less pronounced in the case of $N_{s}=15$, where we still
observe slight oscillations. Additionally, there is a noticeable overall
improvement in the samples compared to the previous case. However, we should
not anticipate a significant difference in uncertainty between nodes $L$, $T$,
and $X$. Interestingly, unlike the previous case, we observe a better
distribution at node $T$ compared to $R$.

Finally, for $N_{s}=20$, we observe that all nodes except $M$ exhibit similar
magnitudes. Once more, there is a notable overall improvement compared to the
previous case, with nodes $L,T,$ and $X$ showing similar levels of
uncertainty, as well as nodes $N$ and $R$.

Intuition for the observations above can be drawn from the figures 2, 3, and
4, where the synthetic data was generated. For $N_{s}=8$ and point $L$ we can
see how in this narrow peak, only the second temporal node contributes for the
component $c_{1}$. It implies that even that the measurement has larger
magnitude than node $N$ or $R$, then the node $L$ could not be able to capture
the information with perturbations in the parameter space. The same phenomena
explains why $X$ and $M$ are the best and the worst design variables
respectively.

As we increase the number of temporal observations, our measurements can
effectively capture information, even in the presence of potential
perturbations in the parameter space. For $N_{s}=15,20$, the dominant factor
is $c^{Feed}$; when $\tau^{inj}$ is fixed, the solution structure in the PDE
remains consistent, leading to increased peak magnitudes.

Figure 2: Solution to the EDM with $N_{s}=8$ over the design points L, M, N,
P, R, T, X and Y.

Figure 3: Solution to the EDM with $N_{s}=15$ over the design points L, M, N,
P, R, T, X and Y.

Figure 4: Solution to the EDM with $N_{s}=20$ over the design points L, M, N,
P, R, T, X and Y.

Figure 5: Utility functions (time injection $\tau^{inj}$ in $x$ axis and
concentration $c^{Feed}$ in $y$ axis). Left: surrogate models for 8, 15, and
20 temporal nodes of the sensor respectively. Middle: the utility function
with the true model. Right: absolute errors.

Figure 6: Samples of MCMC. In the left surrogate models for $N_{s}=8$. In the
right the samples with the true model

.

Figure 7: Samples of MCMC. In the left surrogate models for $N_{s}=15$. In the
right the samples with the true model

.

Figure 8: Samples of MCMC. In the left surrogate models for $N_{s}=20$. In the
right the samples with the true model

.

High-res model | Surrogate model | Training
---|---|---
$5.6\times 10^{-1}$ | $1.25\times 10^{-4}$ | $6.19\times 10^{2}$
$5.6\times 10^{3}$ | $1.25\times 10^{0}$ |
$4.48\times 10^{4}$ | $1\times 10^{1}$ |
Table 3: Average time of evaluations in seconds. Row one: the forward mapping
$G(\theta,d)$, row two: the utility function $U(d)$ at one design point, row
three: MCMC parameter chain of length 80 000.

##  5 Conclusions

Our main objective in this study was to develop feasible computational methods
for the optimal experimental design in the Equilibrium Dispersive Model (EDM)
for chromatography. To achieve this, we implemented a surrogate model and
analyzed the utility function with varying numbers of temporal measurement
times, aiming to improve both the efficiency and accuracy of parameter
estimation.

The surrogate model, based on PSLI, significantly reduced evaluation time
compared to solving the PDE system with the Koren scheme. Moreover, it
reproduced evaluations from the original model to a good approximation, as
confirmed through simulations of the expected utility functions with different
designs. This observations was enforced by considering accuracy in temporal
profile prediction and simulation of posterior ensembles to compare
concentrations. We note that other methods tested, such as polynomial chaos
expansion, failed to produce satisfactory results due to high gradient zones
in the concentration profiles.

The analysis of utility function plots revealed that beyond 15 observation
time instances there was no further improvement in parameter estimation when
varying initial concentrations and injection times, supposing that both design
variables exceed a minimal threshold value.

## Acknowledgement

This work has been supported by the Research Council of Finland (RCoF) through
the Flagship of advanced mathematics for sensing, imaging and modelling,
decision number 358 944. Moreover, TH was supported through RCoF decision
numbers 353 094 and 348 504. JRRG was supported by the LUT Doctoral School.
The authors thank Vesa Kaarnioja for discussions on the PSLI method.

## References

  * [1] A. Alexanderian.  Optimal experimental design for infinite-dimensional bayesian inverse problems governed by pdes: A review.  arXiv preprint arXiv:2005.12998, 2020.
  * [2] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas.  A-optimal design of experiments for infinite-dimensional bayesian linear inverse problems with regularized $\ell_{0}$-sparsification.  SIAM Journal on Scientific Computing, 36(5):A2122–A2148, 2014.
  * [3] A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas.  A fast and scalable method for a-optimal design of experiments for infinite-dimensional bayesian nonlinear inverse problems.  SIAM Journal on Scientific Computing, 38(1):A243–A272, 2016.
  * [4] A. Alexanderian and A. K. Saibaba.  Efficient d-optimal design of experiments for infinite-dimensional bayesian linear inverse problems.  SIAM Journal on Scientific Computing, 40(5):A2956–A2985, 2018.
  * [5] D. Antos and A. Seidel-Morgenstern.  Continuous step gradient elution for preparative separations.  Separation science and technology, 37(7):1469–1487, 2002.
  * [6] R. C. Aster, B. Borchers, and C. H. Thurber.  Parameter estimation and inverse problems.  Elsevier, 2018.
  * [7] A. C. Atkinson and R. D. Tobias.  Optimal experimental design in chromatography.  Journal of Chromatography A, 1177(1):1–11, 2008.
  * [8] V. Barthelmann, E. Novak, and K. Ritter.  High dimensional polynomial interpolation on sparse grids.  Advances in Computational Mathematics, 12(4):273–288, 2000.
  * [9] J. Beck, B. M. Dia, L. F. Espath, Q. Long, and R. Tempone.  Fast bayesian experimental design: Laplace-based importance sampling for the expected information gain.  Computer Methods in Applied Mechanics and Engineering, 334:523–553, 2018.
  * [10] G. E. Box and H. L. Lucas.  Design of experiments in non-linear situations.  Biometrika, 46(1/2):77–90, 1959.
  * [11] F. Bozorgnia and S. Seyed Allaei.  Numerical approximations of chromatographic models.  Numerical Methods for Partial Differential Equations, 2020.
  * [12] T. Briskot, F. Stückler, F. Wittkopp, C. Williams, J. Yang, S. Konrad, K. Doninger, J. Griesbach, M. Bennecke, S. Hepbildikler, et al.  Prediction uncertainty assessment of chromatography models using bayesian inference.  Journal of Chromatography A, 1587:101–110, 2019.
  * [13] H. Bungartz.  Finite Elements of Higher Order on Sparse Grids.  Berichte aus der Informatik. Shaker, 1998.
  * [14] H.-J. Bungartz and M. Griebel.  Sparse grids.  Acta numerica, 13:147–269, 2004.
  * [15] A. G. Carlon, B. M. Dia, L. Espath, R. H. Lopez, and R. Tempone.  Nesterov-aided stochastic gradient methods using laplace approximation for bayesian design optimization.  Computer Methods in Applied Mechanics and Engineering, 363:112909, 2020.
  * [16] K. Chaloner and I. Verdinelli.  Bayesian experimental design: A review.  Statistical science, pages 273–304, 1995.
  * [17] J. A. Christen and C. Fox.  A general purpose sampling algorithm for continuous distributions (the t-walk).  Bayesian Analysis, 5(2):263–281, 2010.
  * [18] D.-L. Duong, T. Helin, and J. R. Rojo-Garcia.  Stability estimates for the expected utility in bayesian optimal experimental design.  Inverse Problems, 39(12):125008, 2023.
  * [19] A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  Bayesian data analysis.  Chapman and Hall/CRC, 1995.
  * [20] G. Guiochon, D. G. Shirazi, and A. Felinger.  Fundamentals of preparative and nonlinear chromatography.  Academic press, 2006.
  * [21] H. Haario, M. Laine, A. Mira, and E. Saksman.  Dram: efficient adaptive mcmc.  Statistics and computing, 16(4):339–354, 2006.
  * [22] H. Haario, E. Saksman, and J. Tamminen.  Adaptive proposal distribution for random walk metropolis algorithm.  Computational Statistics, 14(3):375–395, 1999.
  * [23] H. Haario, E. Saksman, and J. Tamminen.  An adaptive metropolis algorithm.  Bernoulli, pages 223–242, 2001.
  * [24] Q.-L. He and L. Zhao.  Bayesian inference based process design and uncertainty analysis of simulated moving bed chromatographic systems.  Separation and Purification Technology, 246:116856, 2020.
  * [25] T. Helin, N. Hyvonen, and J.-P. Puska.  Edge-promoting adaptive bayesian experimental design for x-ray imaging.  SIAM Journal on Scientific Computing, 44(3):B506–B530, 2022.
  * [26] W. Heymann, J. Glaser, F. Schlegel, W. Johnson, P. Rolandi, and E. von Lieres.  Advanced error modeling and bayesian uncertainty quantification in mechanistic liquid chromatography modeling.  Journal of Chromatography A, 1708:464329, 2023.
  * [27] D. B. Hibbert.  Experimental design in chromatography: A tutorial review.  Journal of Chromatography B, 910:2–13, 2012.  Chemometrics in Chromatography.
  * [28] X. Huan.  Accelerated bayesian experimental design for chemical kinetic models.  Master’s thesis, Massachusetts Institute of Technology, 2010.
  * [29] X. Huan.  Numerical approaches for sequential Bayesian optimal experimental design.  PhD thesis, Massachusetts Institute of Technology, 2015.
  * [30] X. Huan and Y. M. Marzouk.  Simulation-based optimal bayesian experimental design for nonlinear systems.  Journal of Computational Physics, 232(1):288–317, 2013.
  * [31] S. Javeed.  Analysis and Numerical Investigation of Dynamic Models for Liquid Chromatography.  PhD thesis, Otto von Guericke University Magdeburg, 2013.
  * [32] S. Javeed, S. Qamar, W. Ashraf, G. Warnecke, and A. Seidel-Morgenstern.  Analysis and numerical investigation of two dynamic models for liquid chromatography.  Chemical Engineering Science, 90:17–31, 2013.
  * [33] S. Javeed, S. Qamar, A. Seidel-Morgenstern, and G. Warnecke.  Efficient and accurate numerical simulation of nonlinear chromatographic processes.  Computers & chemical engineering, 35(11):2294–2305, 2011.
  * [34] K. Kaczmarski.  Estimation of adsorption isotherm parameters with inverse method—possible problems.  Journal of Chromatography A, 1176(1-2):57–68, 2007.
  * [35] A. Klimke.  Uncertainty Modeling using Fuzzy Arithmetic and Sparse Grids.  PhD thesis, Universität Stuttgart, 2005.
  * [36] A. Klimke.  Efficient construction of hierarchical polynomial sparse grid interpolants using the fast discrete cosine transform.  Inst. for Applied Analysis and Numerical Simulation, Univ. of Stuttgart, Stuttgart, Germany, 7, 2006.
  * [37] A. Klimke.  Sparse grid interpolation toolbox user’s guide, v5. 1.  Inst. for Applied Analysis and Numerical Simulation, Univ. of Stuttgart, Stuttgart, Germany, 2008.
  * [38] A. Klimke and B. Wohlmuth.  Algorithm 847: spinterp: Piecewise multilinear hierarchical sparse grid interpolation in matlab.  ACM Transactions on Mathematical Software (TOMS), 31(4):561–579, 2005.
  * [39] B. Koren.  A robust upwind discretization method for advection, diffusion and source terms.  Centrum voor Wiskunde en Informatica Amsterdam, 1993.
  * [40] Ł. Kubik, R. Kaliszan, and P. Wiczling.  Analysis of isocratic-chromatographic-retention data using bayesian multilevel modeling.  Analytical chemistry, 90(22):13670–13679, 2018.
  * [41] O. A. Ladyzhenskaia, V. A. Solonnikov, and N. N. Ural’tseva.  Linear and Quasi-linear Equations of Parabolic Type, volume 23 of Translations of mathematical monographs.  American Mathematical Society, 1988.
  * [42] O. Le Maître and O. M. Knio.  Spectral methods for uncertainty quantification: with applications to computational fluid dynamics.  Springer Science & Business Media, 2010.
  * [43] M. Leśko, D. Åsberg, M. Enmark, J. Samuelsson, T. Fornstedt, and K. Kaczmarski.  Choice of model for estimation of adsorption isotherm parameters in gradient elution preparative liquid chromatography.  Chromatographia, 78:1293–1297, 2015.
  * [44] G. Lieberman.  Second Order Parabolic Differential Equations.  World Scientific Publishing Company, 2005.
  * [45] D. V. Lindley.  On a measure of the information provided by an experiment.  The Annals of Mathematical Statistics, 27(4):986–1005, 1956.
  * [46] Q. Long, M. Motamed, and R. Tempone.  Fast bayesian optimal experimental design for seismic source inversion.  Computer Methods in Applied Mechanics and Engineering, 291:123–145, 2015.
  * [47] Q. Long, M. Scavino, R. Tempone, and S. Wang.  Fast estimation of expected information gains for bayesian experimental designs based on laplace approximations.  Computer Methods in Applied Mechanics and Engineering, 259:24–39, 2013.
  * [48] Q. Long, M. Scavino, R. Tempone, and S. Wang.  A laplace method for under-determined bayesian optimal experimental designs.  Computer Methods in Applied Mechanics and Engineering, 285:849–876, 2015.
  * [49] J. López-Fidalgo.  Optimal Experimental Design: A Concise Introduction for Researchers, volume 226.  Springer Nature, 2023.
  * [50] A. Nestler and G. Thäter.  Parameter estimation and optimal experimental design for the preparative chromatography.  Numerical Simulation, Optimization and High Performance Computing, page 50, 2011.
  * [51] E. Novak and K. Ritter.  High dimensional integration of smooth functions over cubes.  Numerische Mathematik, 75(1):79–97, 1996.
  * [52] W. Paweł, K. Łukasz, and K. Roman.  Maximum a posteriori bayesian estimation of chromatographic parameters by limited number of experiments.  Analytical chemistry, 87(14):7241–7249, 2015.
  * [53] C. Robert and G. Casella.  Monte Carlo statistical methods, volume 2.  Springer, 2004.
  * [54] K. J. Ryan.  Estimating expected information gains for experimental designs with application to the random fatigue-limit model.  Journal of Computational and Graphical Statistics, 12(3):585–603, 2003.
  * [55] A. Seidel-Morgenstern.  Experimental determination of single solute and competitive adsorption isotherms.  Journal of Chromatography A, 1037(1-2):255–272, 2004.
  * [56] S. A. Smolyak.  Quadrature and interpolation formulas for tensor products of certain classes of functions.  In Doklady Akademii Nauk, volume 148, pages 1042–1045. Russian Academy of Sciences, 1963.
  * [57] A. Solonen, H. Haario, and M. Laine.  Simulation-based optimal design using a response variance criterion.  Journal of Computational and Graphical Statistics, 21(1):234–252, 2012.
  * [58] P. K. Sweby.  High resolution schemes using flux limiters for hyperbolic conservation laws.  SIAM journal on numerical analysis, 21(5):995–1011, 1984.
  * [59] C. R. Vogel.  Computational methods for inverse problems.  SIAM, 2002.
  * [60] G. W. Wasilkowski and H. Wozniakowski.  Explicit cost bounds of algorithms for multivariate tensor product problems.  Journal of Complexity, 11(1):1–56, 1995.
  * [61] P. Wiczling, A. Bartkowska-Śniatkowska, O. Szerkus, D. Siluk, J. Rosada-Kurasińska, J. Warzybok, A. Borsuk, R. Kaliszan, E. Grześkowiak, and A. Bienert.  The pharmacokinetics of dexmedetomidine during long-term infusion in critically ill pediatric patients. a bayesian approach with informative priors.  Journal of pharmacokinetics and pharmacodynamics, 43:315–324, 2016\.
  * [62] P. Wiczling, A. Kamedulska, and Ł. Kubik.  Application of bayesian multilevel modeling in the quantitative structure–retention relationship studies of heterogeneous compounds.  Analytical Chemistry, 93(18):6961–6971, 2021.
  * [63] D. Xiu.  Numerical methods for stochastic computations.  Princeton university press, 2010.
  * [64] Y. Yamamoto, T. Yajima, and Y. Kawajiri.  Uncertainty quantification for chromatography model parameters by bayesian inference using sequential monte carlo method.  Chemical Engineering Research and Design, 175:223–237, 2021.
  * [65] Y. Zhang, G. Lin, M. Gulliksson, P. Forssén, T. Fornstedt, and X. Cheng.  An adjoint method in inverse problems of chromatography.  Inverse Problems in Science and Engineering, 25(8):1112–1137, 2017\.


# The social graph based on real data111Preprint of paper presented at
ICCS’23, https://doi.org/10.1007/978-3-031-36027-5_1

Tomasz M. Gwizdałła<EMAIL_ADDRESS>Faculty of Physics and
Applied Informatics,
University of Lodz, Lodz, Poland  Aleksandra Piecuch  Faculty of Physics and
Applied Informatics,
University of Lodz, Lodz, Poland

###### Abstract

In this paper, we propose a model enabling the creation of a social graph
corresponding to real society. The procedure uses data describing the real
social relations in the community, like marital status or number of kids.
Results show the power-law behavior of the distribution of links and, typical
for small worlds, the independence of the clustering coefficient on the size
of the graph.

##  1 Introduction and Model

The computational study of different social processes requires knowledge of
the structure of the community studied. Among the seminal papers devoted to
study of graph representations of communities we have especially to mention
the Pastor-Satorras paper [10] where different properties of graphs are
studied in the context of epidemic. Now, a lot of problems can be studied with
the use of social graphs. We can mention here the processes which can be
generally called social contagion [1]. Particularly, such topics like disease
spreading [6, 5, 4], transfer of information [9], opinion formation [8, 11] or
innovation emergence [12] are studied with the use of particular graphs. It is
important to emphasize that speaking about social network in the context of
studied problem we have take into account the real community which structure
does not necessarily correspond to the structures obtained by studying,
recently very popular, online networks. The technique which have recently
become popular for such networks is the study of homophily effect [13, 14]. In
our paper, we will present a new approach to creating a society graph, based
on the real statistics describing it. The special attention is devoted to
study the real social links like families, acquaintances. Due to ease of
access, we use the data provided for Polish society, mainly our city - Lodz.

The population is typically presented as an undirected, connected graph, where
nodes correspond to individuals and edges to relations between them. Both
features mentioned above are essential for the understanding of the model.
When assuming the graph’s undirected character, we consider the mutual type of
interaction. It means that compared to well-known social networks, we consider
rather a connection (like the most popular way of using linkedin) than
observation (like twitter). When assuming the connected graph, we ensure that,
in society, we do not have individuals that are completely excluded. We will
come back to this issue when discussing the levels of connections.

We propose to use initially four levels of connection between
individuals/nodes

  * •

I. Household members and inmates. In this class, we collect the people who
have permanent contact - members of families, cohabitants, or people who rent
apartments together.

  * •

II. In this class, several types of connections are considered. The crucial
property of people in these groups is that they are strongly connected but
outside households. We include small groups created in the school or work
environment and further family contacts here. These people meet themselves
every day or almost every day, and small social distances characterize their
meetings.

  * •

III. Acquaitances. In this class, we consider loosely related individuals but
having the possibility of frequent contacts - kids in schools or employees in
companies.

  * •

IV. Accidental but possible contacts - here we sample mainly the people living
in the common areas - villages, city districts. We can meet them accidentally
on the street, in the shop, but these meetings are strongly random.

In the paper, we focus on the first two categories. However, the creation of a
third class will also be described because the second one is based on it. We
assume that at least one individual is always considered as a family member
(1st class) or caregiver (2nd class).

To ensure the correctness of the data, we use only the ones which follow some
basic criteria.

  * •

The source of data can be considered reliable. This criterion, indeed, led to
the strong restriction that all data come from the polish GUS (Statistics
Poland - https://stat.gov.pl)

  * •

The studies were conducted in a similar period. It is tough to find all data
sampled simultaneously during the same study because they cover different
aspects of social or economic statistics.

  * •

The studies were conducted on a similar statistical population. Sometimes, the
easily accessible data concern different communities, e.g., big cities and
village populations. In such a situation, we do rather consider the global
data for the whole country (like the employment vs. size of company
statistics).

In the paper, we model the closed community corresponding rather to a big city
than to scattered areas of villages and small towns. Such a choice is made for
several reasons. The regulations for such environments are more regular,
without many exceptions, like e.g., the size of schools. The modeling of the
distribution of individuals, for example, between districts, is also easier
than creating a large number of small communities, where additional conditions
should be taken into account (like e.g., the general possibilities of
connecting between such communities). The data presented below differ
significantly for various communities. In the calculations, we mainly use the
data for our city - Lodz (Łódź), or for the whole of Poland.

The successive steps of the proposed approach is as follows:

  * •

We start by selecting the correct number of population members according to
the age pyramid. Here we use the data for Lodz, and the borders of intervals
for which data are known are given by
$\\{2,6,12,15,18,24,34,44,54,64,75,95\\}$. For a relatively small percentage
of people in the last interval (above 95), we merge the last two intervals.

  * •

The marital status of men is sampled. We use the data concerning the five
possible states: bachelor, married, cohabitant, divorced, and a widower for
different age intervals (different from the age pyramid). We start from the
men subset because the total number of men is smaller than the number of
women.

  * •

The marriages and cohabitation relations are created as the first level
connection. We use the distribution of differences between men’s and women’s
ages in relations and sample the partners from the women’s subset. For the
remaining women, we attach the status in the same manner as for men. The
correctness of the final marital distribution is the test for the correctness
of the whole procedure, but we do not present the results confirming this
effect here.

  * •

The kids are attached as the first-level connection. Having the distribution
of kids typical for a particular age and the type of family (pair, single men,
single woman), we attach the number of kids of an age different from the age
of the woman (or man) by a value from interval [18,40].

  * •

The kids and youth are attached to schools as the 3rd-level connections. We
know the average size of a school of a particular type, so we can, taking into
account the number of individuals in the corresponding age, determine the
number of schools. Then we create cliques, attaching every individual under 18
(in accordance with Polish law) to the one of the established schools.

  * •

People working in the same companies are connected as the 3rd-level
connections. The basis for assigning people to jobs is the distribution of a
number of employees in one of the four groups specified by polish statistics
(micro, small, average, big - sometimes also the group of large companies is
described, but we omit it). After assigning the type of company to every
individual of the right age (before the retirement age, different for men and
women), we create cliques of size sampled uniformly from the interval
characteristic for particular group.

  * •

The selector of narrower groups as the 2nd-level connections. We choose the
smaller cliques among the 3rd-level connections. We can choose them
arbitrarily, but for the presented calculations, we use the Poisson
distribution with an average equal to $3$ as the size of these cliques.

  * •

Finalization. After the steps described above, there can still be some
unconnected nodes. It mainly concerns retired people, who are not in pairs at
the first level, are not in a household with children and do not work
professionally. For these people, we attach the person, younger by 20-40
years, who can be considered as a separately living child or a caregiver.

The procedure described above allows the creation of the connected graph when
the first and second-level connections are considered. Intentionally, we do
not describe the procedure of creating connections belonging to the 4th-level
because it is not analyzed in the paper.

Figure 1: The comparison of links distribution histograms for different graph
sizes. The solid lines, created for selected set of graph sizes ($10^{3}$,
$10^{4}$,$10^{5}$,$10^{6}$) correspond the power-law dependencies and the
colors correspond to the color of marker for particular size of community.

##  2 Results and Conclusions

We concentrate our attention on four parameters that are typically used to
describe the properties of a graph as the form representing the community.
They are the number of neighbors/links distribution, radius and diameter, and
clustering coefficient. The results are presented for graph sizes from 1000 to
1000000, distributed uniformly in the log scale. The calculations for every
number are performed several times, and the presented values are averaged over
independent runs. The number of repetitions varies from 30 for the smaller
sizes to 5 for greater ones. Due to large computational complexity we do not
present results for radius and diameter for largest graphs. The spread of the
averaged values is presented on plots by the bar corresponding to the standard
deviation value.

The distribution of the number of links is shown in Fig.1. The main
observation is that our procedure does not produce hubs. It is clear since in
the presented scheme, we try to create, as the first and second level
connections, cliques with the size distribution sampled from the Poisson
distribution with a relatively small average - $3$. The crucial property is,
in our opinion, that for our approach, we observe the power-law scaling for
the descending part of the plot. The value of the exponent for the presented
cases shows interesting behavior. They are visibly larger than the typical for
BA value $3$; they start from about $8$ for 1000, seem to exceed $11$, and
start to decrease to about $10$ for 100000. This effect needs, however, more
detailed analysis.

Figure 2: The dependence of radii and diameters on the size of community for
graphs created according to the presented procedure.

The dependence of radii and diameters of graphs on the graph size, shown in
Fig.2, presents the typical, generally logarithmic character. We can observe
that the dispersion of diameters is larger than radii and that the values are
larger than the predicted for scale-free Barabasi-Albert network [3].

Figure 3: The dependence of clustering coefficient on the size of community
for graphs created according to the presented procedure. The values were
obtained in the same wat as those from Fig.2.

A significant effect can be observed in the plot of clustering coefficient 3.
For comparison, on the plot, we also present the curve presented for scale-
free Barabasi-Albert network [7, 2]. ??? It is essential that for our network
we obtain, the expected for small-worlds constancy of the coefficient. Only
very small descend is obseved, but its magnitude is of the order of a few
percent.

The presented results show that the proposed model leads to interesting and
promising effects. The crucial ones are the power-law distribution of links
and the constant value of the clustering coefficient. We can easily indicate
the the directions in which the model should be developed. We have to consider
the appropriate way to add hubs; we have to think about the change of
distribution of close acquaintances (2nd level), and we also think about
another organization of connection levels. However, the changes mentioned
above can be, in our opinion, considered as small corrections and will not
substantially affect the properties shown in the paper.

## References

  * [1] Ferraz de Arruda, G., Petri, G., Rodriguez, P.M., Moreno, Y.: Multistability, intermittency, and hybrid transitions in social contagion models on hypergraphs. Nature Communications 14(1) (2023)
  * [2] Bollobas, B., Riordan, O.: Mathematical results on scale-free random graphs. In: Bormholdt, S., Schuster, A.G. (eds.) Handbook of Graphs and Networks. Wiley (2003)
  * [3] Bollobas, B., Riordan, O.: The diameter of a scale-free random graph. Combinatorica 24, 5–34 (2004)
  * [4] Brattig Correia, R., Barrat, A., Rocha, L.M.: Contact networks have small metric backbones that maintain community structure and are primary transmission subgraphs. PLOS Computational Biology 19(2), 1–23 (02 2023). https://doi.org/10.1371/journal.pcbi.1010854
  * [5] Gwizdałła, T.M.: Viral disease spreading in grouped population. Computer Methods and Programs in Biomedicine 197, 105715 (2020)
  * [6] Karaivanov, A.: A social network model of covid-19. PLOS ONE 15(10), 1–33 (10 2020). https://doi.org/10.1371/journal.pone.0240878
  * [7] Klemm, K., Eguíluz, V.M.: Growing scale-free networks with small-world behavior. Phys. Rev. E 65, 057102 (May 2002)
  * [8] Kozitsin, I.V.: A general framework to link theory and empirics in opinion formation models. Scientific Reports 12, 5543 (2022)
  * [9] Kumar, S., Saini, M., Goel, M., Panda, B.: Modeling information diffusion in online social networks using a modified forest-fire model. Journal of Intelligent Information Systems 56, 355–377 (2021)
  * [10] Pastor-Satorras, R., Castellano, C., Van Mieghem, P., Vespignani, A.: Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979 (Aug 2015)
  * [11] Ramaciotti Morales, P., Cointet, J.P., Muñoz Zolotoochin, G., Fernández Peralta, A., Iñiguez, G., Pournaki, A.: Inferring attitudinal spaces in social networks. Social Network Analysis and Mining 13(1) (2023)
  * [12] Sziklai, B.R., Lengyel, B.: Finding early adopters of innovation in social networks. Social Network Analysis and Mining 13(1) (2023)
  * [13] Talaga, S., Nowak, A.: Homophily as a process generating social networks: Insights from social distance attachment model. JASSS 23(2) (2020)
  * [14] Zhao, T., Hu, J., He, P., Fan, H., Lyu, M., King, I.: Exploiting homophily-based implicit social network to improve recommendation performance. In: International Joint Conference on Neural Networks (IJCNN). pp. 2539 – 2547 (2014)


# Ortho-unit polygons can be guarded with at most
$\lfloor\frac{n-4}{8}\rfloor$ guards

J.M. Díaz-Báñez Departamento de Matemática Aplicada II, Universidad de
Sevilla, Spain. Partially supported by Ministerio de Ciencia e Innovación
CIN/AEI/10.13039/501100011033/ (PID2020-114154RB-I00) and European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement #734922. P. Horn Department of Mathematics, University
of Denver. Supported in part by Simons Collaboration Grant 525039 M.A. Lopez
Department of Computer Science, University of Denver<EMAIL_ADDRESS>Supported in part by a University of Denver John Evans Research Award. N.
Marín Posgrado en Ciencia e Ingeniería de la Computación, Universidad Nacional
Autónoma de México<EMAIL_ADDRESS>ORCID: 0000-0002-0222-3254.
Supported by SEP-CONACyT of Mexico. A. Ramírez-Vigueras Instituto de
Matemáticas, UNAM<EMAIL_ADDRESS>Partially supported by PAPIIT
IN105221, Programa de Apoyo a la Investigación e Innovación Tecnologíca UNAM.
O. Solé-Pi Departamento de Matemáticas, Facultad de Ciencias, UNAM.
<EMAIL_ADDRESS>A. Stevens Department of Computer Science,
University of Denver. J. Urrutia Instituto de Matemáticas, Universidad
Nacional Autónoma de México<EMAIL_ADDRESS>ORCID:
0000-0002-4158-5979. Supported in part by PAPIIT IN102117 Programa de Apoyo a
la Investigación e Innovación Tecnológica, UNAM.

###### Abstract

An orthogonal polygon is called an ortho-unit polygon if its vertices have
integer coordinates, and all of its edges have length one. In this paper we
prove that any ortho-unit polygon with $n\geq 12$ vertices can be guarded with
at most $\lfloor\frac{n-4}{8}\rfloor$ guards.

##  1 Introduction

A (simple) polygon $P$ is a sequence of $n$ vertices $v_{1},\ldots,v_{n}$ and
$n$ edges $v_{1}v_{2},v_{2}v_{3},\ldots,v_{n-1}v_{n},v_{n}v_{1}$ such that no
pair of non consecutive edges share a point. In this paper, and following
standard notation in the Art Gallery literature, we will refer to $P$ as the
region enclosed by the union of the vertices and edges of $P$; the boundary of
$P$, denoted as $\partial P$, consists of the edges and vertices of $P$.
Following O’Rourke’s terminology [14], we will assume that $\partial P\subset
P$. A point $p\in P$ guards a point $q\in P$ if the line segment $pq$ joining
them is contained in $P$. A set $S$ of points in $P$ _guards_ $P$ if any point
in $P$ is guarded by at least one point in $S$. A polygon is called an
_orthogonal polygon_ if its edges are horizontal or vertical line segments.

In this paper we will study guarding problems for _ortho-unit polygons_ , that
is orthogonal polygons all of whose edges have length one. Our main result is
that $\lfloor\frac{n-4}{8}\rfloor$ points are always sufficient and sometimes
necessary to guard any ortho-unit polygon.

##  2 Previous work

The study of _guarding_ or illumination problems of polygons on the plane is a
classic area of research in computational geometry. In 1973 V. Klee posed the
problem of determining the minimum number of guards needed to guard any art
gallery whose boundary is represented by a polygon (Hosenberg 1976 [9]).
Klee’s problem was solved by Chvátal [4], who proved that
$\lfloor\frac{n}{3}\rfloor$ guards are always sufficient and sometimes
necessary. A short and elegant proof was later given by Fisk [6] using
triangulations of polygons.

$\cdots$ (a)

$\cdots$ (b)

Figure 1: (a) is a family of polygons in which $\lfloor\frac{n}{3}\rfloor$
guards are required. (b) is a family of orthogonal polygons in which
$\lfloor\frac{n}{4}\rfloor$ guards are required.

For orthogonal polygons, Kahn, Klawa and Kleitman’s [10] proved that
$\lfloor\frac{n}{4}\rfloor$ points are always sufficient to guard any
orthogonal polygon with $n$ vertices. This bound is also tight. Since then a
plethora of papers studying several variations to Klee’s original question
have appeared; see [14, 15, 18].

Orthogonal polygons have been widely studied, among other reasons,since they
model in a more realistic way traditional buildings, which are better modeled
by orthogonal polygons. To our knowledge there are at least six different
proofs of the Orthogonal Art Gallery Theorem, e.g. see [7, 13, 16, 17]. Hofman
and Kriegel [8] proved that any orthogonal polygon with holes can always be
guarded with $\lfloor\frac{n}{3}\rfloor$ vertex guards. Bjorling-Sachs proved
that $\lfloor\frac{3n+4}{6}\rfloor$ edge guards are always sufficient and
sometimes necessary to guard an orthogonal polygon [3]. Katz and Roisman [12]
studied the problem of guarding the vertices of an orthogonal polygon is
$N\\!P$-hard.

The guarding problem for polyominoes, orthogonal polygons that are obtained by
joining sets of unit squares, was studied by Biedl et al. [2]. They proved
that $\lfloor\frac{m+1}{3}\rfloor$ point guards are always sufficient and
sometimes necessary to guard a polyomino with area $m$, possibly with holes.
Katz and Morgenstern [11] studied the problem of guarding an orthogonal art
gallery with cameras that slide along horizontal or vertical straight tracks.
They proved that if only horizontal tracks are allowed, then finding the
minimum number of tracks needed can be found in polynomial time. Durocher et
al. [5] extended Katz and Morgenstern results by proving that finding a set of
tracks that minimizes the sum of their lengths can also be solved in
polynomial time, even for orthogonal polygons with holes.

Biedl et al. [1] studied the problem of guarding orthogonal polygons with
sliding $k$-transmitters, a point guard that slides along an axis line segment
$s$. A $k$ transmitter can see s point $p\in P$ if the perpendicular line
segment joining it to $S$ intersects the boundary of $P$ at most $k$ times.

The interested reader can find more information in a book by O’Rourke [14] and
two survey papers, by Shermer [15] and Urrutia [18].

##  3 Guarding ortho-unit polygons 1

In this section we prove our main result, namely:

######  Theorem 1.

Any ortho-unit polygon with $n\geq 12$ vertices can be guarded with at most
$\lfloor\frac{n-4}{8}\rfloor$ guards. Moreover, $\lfloor\frac{n-4}{8}\rfloor$
guards is sometimes necessary.

(a)

(b)

Figure 2: (a) is a short ‘macuahuitl’ polygon, with 28 sides, requiring 3
guards. (b) is an illustration of extending a ‘macuahuitl’ polygon. This
requires 8 additional sides and one additional guard. Both polygons require
$\frac{n-4}{8}$ guards.

Suppose $P$ is an ortho-unit polygon with $n\geq 12$ vertices. In order to
prove that $\lfloor\frac{n-4}{8}\rfloor$ guards suffice to guard $P$ we first
introduce an auxiliary structure to help us find appropriate guard positions,
the $\emph{column tree}.$

A _vertical (resp. horizontal) cut_ of $P$ is the extension of a vertical
(resp. horizontal) edge of $P$ at a reflex vertex towards the interior of $P$.
Because $P$ is ortho-unit, each vertical cut connects two reflex vertices. As
observed by O’Rourke in [14], there are $\frac{n-4}{2}$ reflex vertices in any
orthogonal polygon and hence there are $\frac{n}{4}-1$ vertical cuts. They, in
turn, partition $P$ into $\frac{n}{4}$ one wide rectangles, which we call
columns. Note that the top and bottom of each column are horizontal unit
length edges.

abcedfghijkl (a) Partition of the polygon into columns. Columns labeled
$a,c,d,j,k,l$ are called teeth.

abcedfghijkl (b) The corresponding column tree. Vertices representing teeth
are leaves in the column tree.

Figure 3: When two columns share a border in figure (a), their corresponding
vertices share an edge in figure (b).

The _column tree_ of $P$, denoted by $T$, is the tree whose vertices are the
$\frac{n}{4}$ columns of $P$, where two vertices of $T$ are joined by an edge
if they are separated by a vertical cut. Note that a guard placed on a
vertical cut will guard both incident columns. A simple observation, then, is
that a perfect matching in $T$ would yield a set of $\frac{n}{8}$ guards which
guard all of $P$: each of the edges in the matching corresponds to a vertical
cut, and placing a guard on each of these vertical cuts would result in a set
of guards which guard the entire polygon. Of course, $T$ will typically not
contain a perfect matching but this observation is our starting point,
nonetheless.

The following proposition, which follows easily by induction, reveals that
leaves of $T$ are the barrier to finding a perfect matching. We state it in a
form which will be convenient to us later in the proof.

######  Proposition 2.

Suppose $T$ is a tree, and $x\in V(T)$ is a leaf. Then there is a matching in
$T$ which saturates $x$ and all interior vertices of $T$.

Because $P$ is ortho-unit, leaves in $T$ correspond to columns in $P$ that are
one unit high, and have edges to their left or right. We call such structures
_teeth_ in $P$; teeth can be either left- or right-facing depending on whether
they have edges to their right or left side, respectively. As a second
auxiliary structure, we create the _tooth graph_ $H$ of $P$ as follows: the
vertices of $H$ are the teeth of $P$. Two teeth $t_{1}$ and $t_{2}$ are
connected by an edge in $H$ if the horizontal cut from either the top or
bottom edge of $t_{1}$ coincides with a horizontal cut from the top or bottom
edge of $t_{2}$. Equivalently, two teeth are connected by an edge if a single
guard can simultaneously view both.

######  Lemma 3.

The components of $H$ consist of paths and isolated vertices.

###### Proof.

The maximum degree of a tooth in $H$ is clearly two as there are only two
horizontal cuts incident to it, one each coming from its top and bottom edge.
Since the teeth have height one, either _both_ the top and bottom horizontal
cuts of two teeth coincide, or the top cut of one tooth coincides with the
bottom cut of its neighbor. This immediately rules out cycles, and so each
component is a path. ∎

An _interior edge_ of $T$ is an edge of $T$ which is not incident to a leaf.
Each interior edge corresponds to a vertical cut which is _not_ incident to a
tooth. We say that an interior edge of $T$ can guard a tooth if a guard can be
placed on the corresponding vertical cut in such a way that it guards the
tooth.

The following lemma is the key to efficiently guarding the teeth of the
polygon.

######  Lemma 4.

Suppose $P$ is an ortho-unit polygon. Let $O$ denote the set of odd components
of the tooth graph $H$ of $P$. Then there is a matching between $O$ and the
interior edges of $T$ which saturates $O$ and so that each matched interior
edge can guard a leaf in the corresponding matched component.

abcdefg Figure 4: Above is a matching of teeth. We see that $c$ and $g$ are
matched, as are $b$ and $e$ and, $a$ and $d$. $f$ is unmatched.

###### Proof.

Components in $O$ are either isolated teeth or odd length paths of teeth. Note
that the odd length paths can be arranged vertically, and there is a top-most
and bottom-most edge of the top and bottom teeth of the path.

We build a bipartite graph where one partite set consists of components in $O$
and the other partite set consists of interior edges of $T$ as follows. Fix a
component $x$ in $O$. This component has a top-most and bottom-most horizontal
edge (if the component consists of only one tooth, then these are the top and
bottom edge of that tooth, otherwise they lie in the top and bottom tooth
respectively. Consider the horizontal cut defined by the top-most edge; this
intersects $P$ again at a reflex vertex, $r$. Because this is the topmost edge
in a component of $H$, $r$ does not lie on a tooth. Hence, the vertical cut
through $r$ corresponds to an internal edge $e$ of $T$. We connect the
component $x$ to the internal edge $e$, and similarly connect $x$ to the
internal edge obtained in a similar way using the bottom edge.

abcedfghijklmn (a) A matching of internal edges and teeth using Hall’s
Theorem.

abcedfghijklmn (b) A representation of the matching in the column tree. Bold
edges are unmatched.

The internal edges have maximum degree two, their potential neighbors being
determined by the horizontal cuts through the reflex vertices at the top and
bottom. Meanwhile, the odd components have degree _exactly_ two, unless the
top and bottom of the component connect to the same interior edge (and hence
the component is the only neighbor of the edge.) The desired matching, then,
is an easy consequence of Hall’s theorem. That the edge guards a tooth it is
matched to is a consequence of the fact that the horizontal cut coming from
one of the reflex vertices ending the vertical cut corresponds with a
horizontal cut coming from the tooth. ∎

We think of the matching guaranteed by a Lemma 4 as a matching between leaves
of $T$ (teeth) and internal edges so that one of the endpoints of each odd
component of $H$ is matched.

Now we are ready to turn to the meat of the proof of Theorem 1. For our ortho-
unit polygon $P$, consider a matching guaranteed by Lemma 4. As several
matchings are considered in the proof, we call this one the
$\mathcal{L}$-matching. While the $\mathcal{L}$-matching the proof provides
may match teeth with internal edges far from them, it is easier to deal with a
‘closer’ matching. To that end, we consider a matching satisfying the
conclusion of the lemma that minimizes the sum of distances from the teeth
(leaves) to the matched internal edges. Such a matching has the property that
all edges between a tooth and its matched internal edge are also matched (as
otherwise there is a ‘better’ matching.) Fix such a minimal
$\mathcal{L}$-matching; we will work with it for the remainder of the proof.

abcedfghijklmn (a) A matching which minimized the sum of distances from the
teeth to the matched internal edges.

abcedfghijklmn (b) A representation of the matching in the column tree. Edges
shown bold are unmatched.

Place guards on the vertical cuts corresponding to $\mathcal{L}$-matched
internal edges in such a way that they also guard their matched tooth. Let
$\mathcal{F}^{\prime}$ be the forest spanned by $\mathcal{L}$-matched internal
edges in $T$, and $\mathcal{F}$ be the forest obtained by removing all
vertices in $\mathcal{F}^{\prime}$ along with the $\mathcal{L}$-matched teeth.
The vertices in $\mathcal{F}$, then, correspond to columns that still need to
be guarded. Denote by $H^{\prime}$ the graph obtained by removing all
$\mathcal{L}$-matched teeth from the tooth graph $H$. Note that $H^{\prime}$
consists solely of even paths, and hence contains a perfect matching, which we
call the $\mathcal{T}$-matching.

######  Lemma 5.

Let $X$ be a set of internal edges which span a component of
$\mathcal{F^{\prime}}$. Then the guards placed on the corresponding vertical
cuts guard a total of $2|X|+1$ columns of $P$, including teeth.

###### Proof.

Each guard sees both columns incident to the vertical cut it sits on. Since
the internal edges span a tree whose vertices are these columns, this is
precisely $|X|+1$ columns. Furthermore each guard also guards the tooth
$\mathcal{L}$-matched to its internal edge. Since all edges are internal,
these are disjoint from the columns already counted and $2|X|+1$ columns are
guarded in total. ∎

Now we need to guard the remaining columns corresponding to vertices in
$\mathcal{F}$. To this end, we fix a root in $T$. If $\mathcal{F}^{\prime}$ is
non-empty, we choose the root to be a leaf $\mathcal{L}$-matched to an
internal edge in $\mathcal{F}^{\prime}$. If all components in $H$ are even,
then $\mathcal{F}^{\prime}$ is empty; in this case, we choose the root to be
an arbitrary leaf. We now define, for every component of $\mathcal{F}$, its
root to be the vertex in the component closest to the root of $T$.

Let $\ell$ be a leaf of a component of $\mathcal{F}$. We say $\ell$ _blames_ a
component $C$ of $\mathcal{F}^{\prime}$ if $\ell$ is the first vertex outside
of $C$ on the shortest path from $C$ to the root in $T$. From the definition,
one immediately has

######  Lemma 6.

No component $C$ of $\mathcal{F}^{\prime}$ is blamed by more than one leaf
$\ell$.

afgst Figure 7: Suppose the column tree is rooted at $a$. The red, blue,
yellow and green subtreees are the components of $\mathcal{F}^{\prime}$. The
vertices $f,g,s$ and $t$ are in $\mathcal{F}$. While $t$ blames the red
component, $g$ blames both the green and blue components. The yellow component
remains unblamed.

On the other hand,

######  Lemma 7.

Suppose $\ell$ is a leaf or isolated vertex in $\mathcal{F}$ that is not a
leaf in $T$. If $\ell$ is isolated or not the root of its component, then
$\ell$ blames some component $C$ in $\mathcal{F^{\prime}}$.

###### Proof.

Fix a vertex $\ell$ satisfying the hypothesis of the lemma. We claim that
$\ell$ is incident to multiple internal edges of $T$. If $T$ were a star, then
it is easy to see that $P$ is a vertically oriented double comb, and actually
in this case the center vertex is not a leaf in $\mathcal{F}$ (which in the
case is all of $T$, as the components of $H$ in this case are all edges.) Thus
$\ell$ is certainly incident to at least one internal edge of $T$.

If $\ell$ were only incident to one internal edge of $T$ – which is not
$\mathcal{L}$-matched to any tooth, as $\ell$ is in $\mathcal{F}$ and not
$\mathcal{F}^{\prime}$ – then it must also be adjacent to some collection of
leaves. Since $\ell$ is a leaf in $\mathcal{F}$ these leaves must be
$\mathcal{L}$-matched. Then the corresponding tooth could be guarded by the
internal edge incident to $\ell$. Our choice of $\mathcal{L}$-matching then
rules out this case, as we could replace the $\mathcal{L}$-matched internal
edge of some tooth by the internal edge incident to $\ell$ and reduce the
total distance of the $\mathcal{L}$-matching.

Thus $\ell$ is incident to at least two internal edges of $T$. One of these is
not in $\mathcal{F}$ and, since $\ell$ is either isolated in $\mathcal{F}$ or
not the root of its component, not on the shortest path to the root. The other
end of this edge is a vertex which is in $\mathcal{F^{\prime}}$, and $\ell$
blames that vertex’s component as desired. ∎

We are now ready to complete the proof.

###### Proof of Theorem 1.

It remains to guard the vertices of $\mathcal{F}$. For each non-trivial
component, we apply Proposition 2 to obtain a matching in $\mathcal{F}$, the
$\mathcal{F}$-matching, which covers all interior vertices of the component
and the root of each component (taking that root to be $x$ in Proposition 2 if
the root is a leaf). We place guards on each of the vertical cuts
corresponding to these matched edges. At this point, we have placed guards on
vertical cuts the edges of $\mathcal{F^{\prime}}$ and the
$\mathcal{F}$-matched edges of $\mathcal{F}$. The only unobserved columns now
correspond to isolated vertices and unmatched leaves in $\mathcal{F}$.

For each unmatched leaf or isolated vertex in $\mathcal{F}$ that is not a leaf
in $T$ we place a guard dedicated to guarding the associated column.

The remaining leaves correspond to teeth in $P$ and are handled as follows:
fix a perfect $\mathcal{T}$-matching in $H^{\prime}$, the tooth graph
consisting of teeth that were not matched in the first step. For each edge of
this $\mathcal{T}$-matching, we add guards as follows:

  1. 1.

If both teeth are saturated in the $\mathcal{F}$-matching we do nothing, as
both teeth are already guarded.

  2. 2.

If one tooth is saturated in the $\mathcal{F}$-matching but one is not, note
that the guard can be placed in such a way that it actually guards _both_
teeth simultaneously, while still guarding the columns it was originally
responsible for.

  3. 3.

If neither tooth is saturated in the $\mathcal{F}$-matching, place a single
guard responsible for both teeth.

By this construction all columns are guarded. Furthermore, on average each
guard is responsible for at least two unique columns. This follows as the only
guards initially assigned a single column are those which, by Lemma 7 blame
some component of $\mathcal{F}^{\prime}$. Per Lemma 5, the guards guarding the
columns of a component of $\mathcal{F}^{\prime}$ guard, in the aggregate, an
extra column beyond the claimed average of two each. Assigning this ‘extra’
column to the guard only assigned a single column which blames that component
gives the claimed average.

This suffices to prove that at most $n/8$ guards can guard $P$. To prove that,
actually, at most $(n-4)/8$ guards suffice, we proceed as follows.

If $\mathcal{F}^{\prime}$ is non-empty, note that the root of $T$ is incident
to an internal vertex $x$ of $T$ that by our choice of matching is in
$\mathcal{F}^{\prime}$. The component of $x$ is not blamed by any leaf and
hence the ‘extra’ column is never assigned. Thus, if $g$ guards are placed
$2g+1$ is a lower bound for the number of distinct columns guarded so that
$2g+1\leq\frac{n}{4}$; rearranging $g\leq\frac{n-4}{8}$ as desired.

If $\mathcal{F}^{\prime}$ is empty, which occurs if all components of $H$ are
even, there is an ‘extra’ column assigned to a guard if, when considering the
$\mathcal{T}$-matching, case (b) above ever occurs.

If case (b) never occurs, take an alternating path in the matching on
$\mathcal{F}$ starting from some leaf $\ell$ in $T$. This ends at some other
leaf of $T$, $\ell^{\prime}$.

If $\ell$ and $\ell^{\prime}$ are not paired in the $\mathcal{T}$-matching
with each other, interchange the edges of the alternating path to obtain a new
matching. For this matching it is easy to see that case (b) will occur and the
matching satisfies all other desired conditions.

If $\ell$ and $\ell^{\prime}$ are paired in the $\mathcal{T}$-matching, ensure
that we are in case (c) above, possibly interchanging the edges of the
alternating path to obtain a new $\mathcal{F}$-matching. Note that the path
between $\ell$ and $\ell^{\prime}$ in $T$ then contains a matched internal
edge, and a guard can be placed on this internal edge in such a way that it
sees not only the incident columns but also guards _both_ the teeth
corresponding to both $\ell$ and $\ell^{\prime}$.

This assignment of guards to columns ensures that there is at least one more
than twice the number of guards columns guarded, and hence yields the desired
bound.

∎

##  4 Integral Polygons

An orthogonal polygon is an _integral orthogonal polygon_ if all of its edges
have integer length. The lower bound on the number of guards needed to guard
this family of polygons is given by the integral version of the well-known
orthogonal comb polygon which requires $n/4$ guards. See Figure 8 for an
example.

Figure 8: An integral polygon $P$ of perimeter $N$ which requires $N/6$
guards. The elements of a guard set for $P$ are shown in red.

We can prove, via an unenlightening and complicated proof (not included in
this paper), that integral orthogonal polygons of perimeter $N$ can be guarded
using at most $\lfloor N/5\rfloor$ guards. Nevertheless, we conjecture that
$N/6$ guards are always sufficient:

######  Conjecture 8.

Let $P$ be an integral orthogonal polygon of perimeter $N$. Then $P$ can
always be guarded with at most $\lfloor N/6\rfloor$ guards.

## References

  * [1] T. Biedl, T. M. Chan, S. Lee, S. Mehrabi, F. Montecchiani, H. Vosoughpour, and Z. Yu.  Guarding orthogonal art galleries with sliding k-transmitters: Hardness and approximation.  Algorithmica, 81(1):69–97, 2019.
  * [2] T. Biedl, M. T. Irfan, J. Iwerks, J. Kim, and J. S. Mitchell.  The art gallery theorem for polyominoes.  Discrete & Computational Geometry, 48(3):711–720, 2012.
  * [3] I. Bjorling-Sachs.  Edge guards in rectilinear polygons.  Computational Geometry, 11(2):111–123, 1998.
  * [4] V. Chvátal.  A combinatorial theorem in plane geometry.  Journal of Combinatorial Theory, Series B, 18(1):39–41, 1975.
  * [5] S. Durocher, O. Filtser, R. Fraser, A. D. Mehrabi, and S. Mehrabi.  Guarding orthogonal art galleries with sliding cameras.  Computational Geometry, 65:12–26, 2017.
  * [6] S. Fisk.  A short proof of Chvátal’s watchman theorem.  Journal of Combinatorial Theory, Series B, 24(3):374, 1978.
  * [7] E. Györi.  A short proof of the rectilinear art gallery theorem.  SIAM Journal on Algebraic Discrete Methods, 7(3):452–454, 1986\.
  * [8] F. Hoffmann and K. Kriegel.  A graph-coloring result and its consequences for polygon-guarding problems.  SIAM Journal on Discrete Mathematics, 9(2):210–224, 1996.
  * [9] R. Honsberger.  Mathematical gems ii, mathematical association of america, dolciani math, 1976.
  * [10] J. Kahn, M. Klawe, and D. Kleitman.  Traditional galleries require fewer watchmen.  SIAM Journal on Algebraic Discrete Methods, 4(2):194–206, 1983\.
  * [11] M. J. Katz and G. Morgenstern.  Guarding orthogonal art galleries with sliding cameras.  International Journal of Computational Geometry & Applications, 21(02):241–250, 2011.
  * [12] M. J. Katz and G. S. Roisman.  On guarding the vertices of rectilinear domains.  Computational Geometry, 39(3):219–228, 2008.
  * [13] J. O’Rourke.  An alternate proof of the rectilinear art gallery theorem.  Journal of Geometry, 21(1):118–130, 1983.
  * [14] J. O’Rourke.  Art Gallery Theorems and Algorithms, volume 57.  Oxford University Press, Inc., New York, NY, USA, 1987.
  * [15] T. C. Shermer.  Recent results in art galleries (geometry).  Proceedings of the IEEE, 80(9):1384–1399, 1992.
  * [16] G. T. Toussaint.  Computational morphology: A Computational geometric approach to the analysis of form.  Elsevier, 2014.
  * [17] J. Urrutia.  Sixth proof of the orthogonal art gallery theorem.  Technical Report TR-97-03, Department of Computer Science, University of Ottawa, 1997.
  * [18] J. Urrutia.  Art gallery and illumination problems.  In J.-R. S. Urrutia, editor, Handbook of Computational Geometry, pages 973 – 1027. North-Holland, Amsterdam, 2000.