\section{Introduction}
It has been said that seismology is the most successful geophysical technique for learning about the Earth's interior~\citep[Section~9.4]{AkiRichards2002}.
A common approach used by seismologists to gain information about the stratified subsurface of the Earth is traveltime inversion.

To perform such an inversion, seismologists require at least two mathematical formalisms.
The first is  the forward problem, which, in the case of this paper, is derived using ray theory.
As such, the path of the ray that connects a source and receiver must obey Fermat's principle of stationary time.
The model must account for the material properties of the subsurface, which is satisfied by considering a velocity profile in the form of interval velocities.
The second is the quantitative modelling, which, through a numerical procedure, alters the model parameters of the forward problem to resemble a set of physical observations.
We refer to the second formalism, as the inverse problem, {\it sensu lato}.
The choice of numerical procedure depends greatly on the structure of the forward-problem model.
In this paper, we consider an inverse problem of the least-squares type, which, in the case of traveltime inversion, seeks to minimize the difference between modelled and measured traveltimes.
The optimization technique that we consider is Newton's method in optimization, which is discussed in detail in textbooks of numerical optimization~(e.g.,~\citet[Chapter~6]{NocedalWright1999} and~\citet[Chapter~3]{Fletcher2000}).

In seismology, the combination of these mathematical formalisms forms the basis of many studies.
From~\citet{Dix1955} calculating interval velocities to~\citet{Stewart1984} using least-squares optimization techniques of traveltimes to refine interval-velocity accuracy, seismology has been used routinely to estimate layer properties of the stratified subsurface of the Earth.
Further, \citet{PujolEtAl1985} formulate multilayer ray-tracing equations for traveltime calculation of surface shot gathers.
Therein, they consider a linearized solution for the Taylor series expansion of the expressions of traveltimes and offsets, in terms of the speed and ray parameter, are truncated after the first order.
The resulting linear system is solved using the Levenberg-Marquardt algorithm, which is an iterative approach that uses the information from the gradient to determine the next iterate.
Similarly, \citet{MaoStuart1997} consider a ray-tracing multilayer model with laterally heterogeneous layers for the purpose of transmission and reflection tomography.
Therein, to minimize the misfit of their traveltime expressions, they truncate the Taylor series expansion after the first order.
Consequently, they use the Gauss-Newton algorithm, which approximates the Hessian using first-order derivatives, to minimize the traveltime misfit.

Another specification of ray theory is to model media whose velocity changes not only with depth, but also with direction, which we refer to as an anisotropic inhomogeneous medium.
Indeed, both the anisotropy and inhomogeneity of physical material have an intertwined effect on the propagation of seismic waves.
However, within the framework of mathematics, it is entirely possible to quantify their effects separately.
To this end, there is a wealth of research on the effects of anisotropy and inhomogeneity, of which we highlight~\citet{Cerveny1989,Cerveny2001},~\citet{Winterstein1990},~\citet{Helbig1994},~\citet{Slawinski2020a,Slawinski2020b}.
For the purposes of inversion, it is hopeless to attempt to invert for the model parameters that govern the symmetry class of general anisotropy.
Instead, as indicated by~\citet[p.~91]{Cerveny1989}, it is advisable to consider simpler forms of anisotropy, such as elliptically anisotropic media.

One model of this type is that of~\citet{SlawinskiEtAl2004}, wherein the defining characteristics are that the velocity of a propagating wave increases linearly with depth and exhibits an elliptical velocity dependence along the direction of propagation.
They demonstrate that the statistical significance of the anisotropy parameter, $\chi$\,, is small, but significant, and remain consistent with the weakly anisotropic assumptions discussed by~\citet{Thomsen1986}.
The model is extended by~\citet{RogisterSlawinski2005}, who use Hamilton's equations to obtain the ray-tracing equations and derive analytic solutions for the traveltime equations, which are valid for both upgoing and downgoing ray paths.
This traveltime model has been the topic of subsequent studies during the past decade.

Let us call attention to a few such studies.
\citet{DanekSlawinski2012} use the model to perform a stochastic interpretation of geophysical data to justify the inclusion of the anisotropy parameter by means of the Bayesian Information Criterion.
Therein, to perform the traveltime inversion, they use the Metropolis-Hastings algorithm, which is a random walk algorithm that perturbs model parameters numerically to reduce the traveltime misfit.
In applying a one-layer model to synthetic and real traveltime data, they demonstrate the necessity of the anisotropy parameter to model adequately the observed traveltimes.
\citet{ZhangEtAl2017} approximate the model to calculate numerically ray paths through laterally homogeneous and heterogeneous velocity models.
By imposing triangular grids on the model subsurface, they solve the equations of motion intrinsic to the problem using the Runge-Kutta method.
\citet{GierlachDanek2018} construct two synthetic, multilayer setups, each containing one anisotropic layer.
Using the Nelder-Mead simplex algorithm, they use the traveltime model to invert for the layer properties and, using Bayesian Information Criterion, they determine the optimal number of model parameters to model adequately the traveltimes under the influence of noise.

In this paper, we consider the mathematical model developed by~\citet{SlawinskiEtAl2004}, and extended by~\citet{RogisterSlawinski2005}, to estimate anisotropy and inhomogeneity parameters of the subsurface using vertical seismic profile traveltimes.
Its structure is as follows.
We begin by adapting the traveltime equations for multilayered models, providing an overview of Newton's method and its convergence-ensuring specifications, and establishing the requirement of logarithmic barrier functions to guarantee model-parameter restrictions.
Then, we construct synthetic traveltime models to assess the reliability of results of the optimization under the influence of six noise profiles.
Next, we repeat the optimization for a one- and two-layer model and determine the threshold of noise that yields acceptable results for each model, which we assess using relative errors between the estimated and true model parameters.
Finally, we perform the optimization on a real-data example and conclude with a discussion of the results.
Also, this paper includes an appendix that exemplifies the control experiment of the optimization technique as well as the details of the convergence specifications and parameter restrictions.
\section{Formulation}
\subsection{Mathematical model}
\label{eq:abchiModel}
Let us consider a model of $N$ horizontal layers each exhibiting anisotropy and inhomogeneity.%
\footnote{
	Throughout this paper, in keeping with SI units, the units for $a$ and $b$ are m\,s$^{-1}$ and s$^{-1}$\,, respectively, for traveltime are s\,, for the ray parameter are m$^{-1}$\,s\,
}
The inhomogeneity is a linear increase of velocity, $V(z)=a+b\,z$, where $a$ is the vertical speed of the wavefront at the top of a layer, $b$ is the measure of inhomogeneity within that layer, and $z$ is the vertical component that corresponds to depth.
The anisotropy is an elliptical velocity dependence of a wavefront~\citep{SlawinskiEtAl2004},
\begin{equation}
	\label{eq:chi}
	\chi = \frac{v_{h}^{2}-v_{v}^{2}}{2\,v_{v}^{2}}\,,
\end{equation}
where $v_h$ and $v_v$ are the horizontal and vertical velocity; for $v_h = v_v$, $\chi = 0$ and, hence, the wavefront velocity is isotropic.

We consider a series of $M$ point sources aligned in the $xz$-plane along the surface of the multilayered model and a single receiver.
Along the $j$th ray, the coordinates at the $\ell$th layer interface are given by the sequence
\begin{equation}
S_j = \{(x_{\ell,j},z_\ell)\}_{\ell=1}^{N+1}
\quad{\rm for}\quad j\in\{1\,,\dots,M\}
\,;
\end{equation}
we denote $\boldsymbol{S}$ as the $M\times1$ vector of sequences.
As such, the source coordinates are $(x_{1,j},z_1)$\,, where $x_{1,j}$ are the horizontal offsets and the vertical position $z_1 = 0$ for all $M$ sources.
The receiver coordinates are $(x_{N+1,j},z_{N+1})$\,, where $x_{N+1,j} = 0$ by design and $z_{N+1}$ is the depth of the receiver.
For the $i$th layer, $a_{i}$ is the vertical speed at the top of the layer, $b_{i}$ is the linear increase in vertical speed within the layer, and $\chi_{i}$ is the anisotropy of the layer.
Herein, $b_i>0$ corresponds to the increasing compaction of the subsurface with depth, which increases the speed of a propagating wave.
In view of physical considerations of applying the model to shale, we require $\chi_i>0$ due to the increase in the horizontal component of speed of a propagating wave, which is a directional property of shale.

Thus, the traveltime of a wave propagating along the $j$th ray is
\begin{equation}
	\label{eq:t}
	t(S_j;\boldsymbol{\beta})
	=
	\sum\limits_{i=1}^{N}\frac{\arctanh\left(\mathfrak{p}_{i,j}\,b_i\,(x_{i+1,j}-x_{i,j})-\sqrt{1-{\mathfrak{p}_{i,j}}^{2}{a_i}^2\left(1+2\,\chi_i\right)}\right) + \arctanh\sqrt{1-{\mathfrak{p}_{i,j}}^{2}{a_i}^2\left(1+2\,\chi_i\right)}}{b_i}
	\,,
\end{equation}
where $\boldsymbol{\beta} = \left[\left\{a_i\right\}_{i=1}^{N},\left\{b_i\right\}_{i=1}^{N},\left\{\chi_i\right\}_{i=1}^{N}\right]^t$ is the $3N\times1$ vector of model parameters and
\begin{equation}
	\label{eq:p_{i,j}}
	\mathfrak{p}_{i,j}
	=
	\frac{
		2\,(x_{i+1,j}-x_{i,j})
	}{
		\sqrt{
			\left(\left(x_{i+1,j}-x_{i,j}\right)^2+\left(1+2\,\chi_i\right)\left(z_{i+1}-z_i\right)^2\right)
			\left(\left(2\,a_i+b_i\left(z_{i+1}-z_i\right)\right)^2\left(1+2\,\chi_i\right)+{b_i}^2\left(x_{i+1,j}-x_{i,j}\right)^2\right)
		}
	}
	\,,
\end{equation}
is the ray parameter in the $i$th layer, which is a conserved quantity along the ray.
Herein, we denote $\boldsymbol{t}:=\boldsymbol{t}(\boldsymbol{S};\boldsymbol{\beta})$ as the $M\times1$ vector of modelled traveltimes.
Also, let us remark that the symbol $\boldsymbol{\beta}$ should not be confused with $S$-wave speeds, which is a common use of the symbol in seismology.
To calculate the traveltime along the $j$th ray, we require a complete sequence $S_j$\,, for which the coordinates of the $N-1$ intermediate $x_{\ell,j}$ are yet unknown.
As such, for a $N$-layer model, we form the system of equations
\begin{equation}
	\label{eq:p_{i,j}system}
	\left\{
	\begin{array}{rcl}
		p_{1,j}-p_{2,j} &=& 0 \\
		p_{2,j}-p_{3,j} &=& 0 \\
		&\vdots& \\
		p_{N-1,j}-p_{N,j} &=& 0 \\
	\end{array}
	\right.
\end{equation}
and solve numerically for each $x_{\ell,j}$ for $\ell\in\{2,\dots,N\}$\,.
\subsection{Optimization technique}
\label{sec:OptTech}
To perform the optimization, we require a set of $M$ measured traveltimes that we denote by $\boldsymbol{T}$\,, which is a  $M\times1$ vector of traveltimes.
To minimize the residual between measured and modelled traveltimes, we consider the problem
\begin{equation}
	\label{eq:minf}
	\min_{\boldsymbol{\beta}} f(\boldsymbol{\beta})
	\,,
	\quad{\rm where}\quad
	f(\boldsymbol{\beta})
	=
	\sum\limits_{j=1}^{M}\left(T_j - t(S_j\,;\boldsymbol{\beta})\right)^2
\end{equation}
is a twice-differentiable function, and hence, the gradient, $\boldsymbol{g}(\boldsymbol{\beta})$\,, and Hessian, $\boldsymbol{H}(\boldsymbol{\beta})$\,, of $f(\boldsymbol{\beta})$\,, exist and are continuous, whose components are
\begin{equation}
	\label{eq:g_i,H_ij}
	g(\boldsymbol{\beta})_i = \frac{\partial f(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}_i}
	\quad{\rm and}\quad
	H(\boldsymbol{\beta})_{i,j} = \frac{\partial^2 f(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}_i\partial\boldsymbol{\beta}_j}\,;
\end{equation}
the former is a $3N\times1$ vector and the latter is a $3N\times3N$ matrix.
We use Newton's method in optimization, which relies on second-order Taylor series approximation of $f(\boldsymbol{\beta})$ at $\boldsymbol{\beta} = \boldsymbol{\beta}^{(k)}$\,, to form a quadratic approximation to problem~\eqref{eq:minf}, such that
\begin{equation}
	\label{eq:MultivariateVectorTaylor}
	f(\boldsymbol{\beta}^{(k)} + \alpha\,\Delta\boldsymbol{\beta}^{(k)})
	=
	\left.f(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(k)}}
	+
	(\alpha\,\Delta\boldsymbol{\beta}^{(k)})^t\left.\boldsymbol{g}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(k)}}
	+
	\frac{1}{2}(\alpha\,\Delta\boldsymbol{\beta}^{(k)})^t\left.\boldsymbol{H}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(k)}}(\alpha\,\Delta\boldsymbol{\beta}^{(k)})
\,,
\end{equation}
where
\begin{equation}
	\label{eq:DeltaBeta^{(k)}}
	\Delta\boldsymbol{\beta}^{(k)}
	=
	-\left(\left.\boldsymbol{H}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(k)}}\right)^{-1}\left.\boldsymbol{g}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}^{(k)}}
\end{equation}
is the Newton direction, $\alpha\in[0,1]$ is the step size, and $k$ is an iterative superscript.
To obtain the Newton direction, we set the derivative, with respect to $\Delta\boldsymbol{\beta}^{(k)}$\,, of expression~\eqref{eq:MultivariateVectorTaylor} to zero and solve for $\Delta\boldsymbol{\beta}^{(k)}$\,, where we assume $\boldsymbol{H}(\boldsymbol{\beta})$ is positive definite.
Henceforth, we adopt the notation
\begin{equation}
	\label{eq:f^k_g^k_H^k}
	f^{(k)} := \left.f(\boldsymbol{\beta})\right|_{\boldsymbol{\beta} = \boldsymbol{\beta}^{(k)}}
	\,,\quad
	\boldsymbol{g}^{(k)} := \left.\boldsymbol{g}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta} = \boldsymbol{\beta}^{(k)}}
	\,,\quad
	\boldsymbol{H}^{(k)} := \left.\boldsymbol{H}(\boldsymbol{\beta})\right|_{\boldsymbol{\beta} = \boldsymbol{\beta}^{(k)}}
	\,.
\end{equation}

To determine a local solution of problem~\eqref{eq:minf}, we require an algorithm to generate a sequence of iterates~$\{\boldsymbol{\beta}\}_{k=0}^{\infty}$ that minimize the traveltime residual.
The value of~$f$ at each iterate must be such that~$f^{(k+1)} < f^{(k)}$\,, or, in other words, the sequence of iterates must lead to a strict decrease in $f$\,.
The algorithm is terminated when either a sufficient residual value is achieved or the condition of strict decrease is not met.
To generate the sequence of iterates, we consider a line search Algorithm~\ref{algo:LineSearch}, such that
\begin{equation}
\boldsymbol{\beta}^{(k+1)} = \boldsymbol{\beta}^{(k)} + \alpha^{(k)}\Delta\boldsymbol{\beta}^{(k)}
\,,
\end{equation}
which computes a step size $\alpha^{(k)}$ along the direction $\Delta\boldsymbol{\beta}^{(k)}$ to minimize $f$\,.

\begin{algorithm}[h]
	\small
	\SetAlgoLined
	For $\boldsymbol{\beta}^{(k)}$\,, calculate $f^{(k)}$\,, $\boldsymbol{g}^{(k)}$\,, $\boldsymbol{H}^{(k)}$\,, $\Delta\boldsymbol{\beta}^{(k)}$\,.
	Then, set $r=f^{(k)}$\,, $r' = r$\,, $\alpha' = 0$\,, $\Delta\alpha=0.1$\;
	\While{$r' == r$}{
		$\boldsymbol{\beta}' = \boldsymbol{\beta}^{(k)} + \alpha'\,\Delta\boldsymbol{\beta}^{(k)}$\;
		$r = \left.f(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}'}$\;
		\eIf{$r < r'$}{
			$r' = r$\;
		}{
			\eIf{$\alpha' - 2\Delta\alpha > 0$}{
				$\alpha' = \alpha' - 2\Delta\alpha$\;
			}{
				$\alpha' = 0$\;
			}
			$\Delta\alpha = \Delta\alpha/10$\;
			$\boldsymbol{\beta}' = \boldsymbol{\beta}^{(k)} + \alpha'\Delta\boldsymbol{\beta}^{(k)}$\;
			$r = \left.f(\boldsymbol{\beta})\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}'}$\;
			$r'=r$\;
		}
		$\alpha' = \alpha' + \Delta\alpha$\;
		\If{$\Delta\alpha < 1\times10^{-12}$}{
			break\;
		}
		\If{$\alpha' \geq 1$}{
			break\;
		}
	}
	$\alpha^{(k)} = \alpha'$\;
	$\boldsymbol{\beta}^{(k+1)} = \boldsymbol{\beta}^{(k)} + \alpha^{(k)}\Delta\boldsymbol{\beta}^{(k)}$\;
	\caption{Backtracking line search}
	\label{algo:LineSearch}
\end{algorithm}

For Algorithm~\ref{algo:LineSearch} to produce a sequence of iterates that are globally convergent, the search direction $\Delta\boldsymbol{\beta}^{(k)}$ must be a descent direction, which is true if $\boldsymbol{H}^{(k)}$ is positive definite.
However, as indicated by~\citet[p.~136]{NocedalWright1999}, if $\boldsymbol{H}^{(k)}$ is not positive definite, or near singular, the Newton direction may lead to ascending values of $f$ or may be excessively long.
To ensure positive definiteness of $\boldsymbol{H}^{(k)}$\,, we execute the modified Cholesky factorization algorithm on $\boldsymbol{H}^{(k)}$ prior to calculating the Newton direction.
Details of this algorithm are outlined in~\citet[p.~145--149]{NocedalWright1999}, along with its use in this paper in Appendix~\ref{app:Opt_N=1_NoMods} and its necessity in Appendix~\ref{app:Barrier}.

\subsection{Logarithmic barrier function}
\label{sec:LogBar}
The solution to problem~\eqref{eq:minf} is vector $\widehat{\boldsymbol{\beta}}$\,, which is referred to as the estimator of problem~\eqref{eq:minf}.
In unconstrained optimization, there are no restrictions on the model-parameter values that comprise $\widehat{\boldsymbol{\beta}}$\,.
However, to maintain physically meaningful results in fit-for-purpose applications of minimizing traveltime residuals, it is reasonable to impose restrictions on some, if not all, model parameters.
To do so, we introduce a logarithmic barrier function to $f$\,, which changes problem~\eqref{eq:minf} to an inequality-constrained optimization problem.
Thus, we have
\begin{equation}
	\label{eq:minf_c}
	\min_{\boldsymbol{\beta}} f(\boldsymbol{\beta})
	\quad\text{subject to}\quad
	c_i\geq0\,,\quad
	i\in\mathcal{I}\,,
\end{equation}
where $c_i$ is the set of inequality constraints and $\mathcal{I}$ is the set of integers\,.
The objective function for problem~\eqref{eq:minf_c} is
\begin{equation}
	\label{eq:P(beta)}
	P(\boldsymbol{\beta}) = f(\boldsymbol{\beta}) - \sum\limits_{i\in\mathcal{I}}\,\log l_i\,,
\end{equation}
where the summand of the latter term in $P$ is the logarithmic barrier function.
In a similar manner to notation~\eqref{eq:f^k_g^k_H^k}, we let $P^{(k)}$ denote value of the objective function at the $k$th iteration of the optimization---it is analogous to the $k$th value of $f^{(k)}$ for problem~\eqref{eq:minf}.
Furthermore, in the context of problem~\eqref{eq:minf_c}, we abuse notation~\eqref{eq:g_i,H_ij} by stating equivalently that the components of the gradient and Hessian of $P(\boldsymbol{\beta})$ are
\begin{equation}
	g(\boldsymbol{\beta})_i = \frac{\partial P(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}_i}
	\quad{\rm and}\quad
	H(\boldsymbol{\beta})_{i,j} = \frac{\partial^2 P(\boldsymbol{\beta})}{\partial\boldsymbol{\beta}_i\partial\boldsymbol{\beta}_j}
	\,.
\end{equation}

According to~\citet[Section 17.2]{NocedalWright1999}, the barrier functions have the properties such that they are infinite everywhere except in a particular domain, smooth inside the domain, and their value approaches infinity at the boundary of the domain.
To satisfy these properties, we consider the logistic function,
\begin{equation}
	\label{eq:l}
	l(x) = \frac{1}{1+\exp\left(-r\left(x-x_0\right)\right)}
	\,,
\end{equation}
where $r$ is the growth rate and $x_0$ is its center.
For large $r$\,, equation~\eqref{eq:l} resembles a step function about $x_0$ with a minimum value of zero and a maximum value of one.
Since $\log(0) = -\infty$ and $\log(1) = 0$\,, we exploit~\eqref{eq:l} to restrict the values of model parameters.
We provide further details on the use of the logistic function  as a logarithmic barrier function in Appendix~\ref{app:Barrier}.

\section{Synthetic-data testing}
Let us assess the reliability of results in using Newton's method to solve problem~\eqref{eq:P(beta)}.
For consistency, we retain the same source-receiver combinations across all executions of the algorithm.
We specify $M=139$ point sources along the surface, with offsets ranging from approximately 80\,m to 3300\,m\,, at intervals of approximately 25\,m\,, and a single receiver at a depth of approximately 1850\,m\,.
To ensure physically meaningful results, we use the logarithmic barrier function to restrict model parameters such that $b_i>0$ and $\chi_i>0$\,.

First, we set up the following control experiment.
For a known set of model parameters, $\boldsymbol{\beta}$\,, we generate sequences $\boldsymbol{S}$ and calculate traveltimes $\boldsymbol{t}(\boldsymbol{S};\boldsymbol{\beta})$\,, which we designate as the measured traveltimes $\boldsymbol{T}$\,.
Then, from a set of starting values, $\boldsymbol{\beta}^{(1)}$\,, we perform Newton's method to obtain the sequence of iterates $\{\boldsymbol{\beta}^{(k)}\}$ that minimize $f$ and, hence, solve problem~\eqref{eq:P(beta)}, whose solution is the estimator vector, $\widehat{\boldsymbol{\beta}}$\,.
We quantify the ``goodness of fit'' by calculating the relative errors, $\delta(\widehat{\boldsymbol{\beta}}) := ((\widehat{\boldsymbol{\beta}} - \boldsymbol{\beta})/\boldsymbol{\beta})\times100\%$
between the observed $\widehat{\boldsymbol{\beta}}$ and expected $\boldsymbol{\beta}$\,.

Second, we assess the reliability of results under the influence of random noise, which serves to simulate errors in the traveltime data.
The range of the random noise is set at a constant percentage of the measured traveltimes used in the control experiment.
In other words, as the traveltime increases with offset, the range of the noise increases by a fixed percentage.
Within the fixed-percentage range, we generate the noise profile using Matlab's {\tt rand} function, which returns uniformly distributed random numbers within a determined interval.
We generate six noise profiles for which the fixed-percentage range increases by a factor of 10 from $\pm0.000001\%$ to $\pm0.1\%$\,.
As the range increases for each noise profile, the magnitude of noise within the noise profile also increases.
The corresponding mode of the orders of magnitude increases by a factor of 10 from nanoseconds $(10^{-9})$ to milliseconds $(10^{-3})$\,.
Throughout the assessment, we use the same six noise profiles that are plotted in Figure~\ref{fig:NoiseProfiles}.
Similarly to the control experiment, we quantify the goodness of fit of estimator~$\widehat{\boldsymbol{\beta}}$\,.

\begin{figure}[h]
	\centering
	\includegraphics[width=\textwidth]{Fig_NoiseProfiles.pdf}
	\caption{\small%
		Six noise profiles for increasing percentage ranges of traveltimes
	}
	\label{fig:NoiseProfiles}
\end{figure}

Third, we perform the control and noise-assessment experiments for a one- and two-layer model.
This process is repeated on eight sets of model parameters, with increasing anisotropy.
We distinguish one set of model parameters from another using notation $[\boldsymbol{\beta}]_s$ for $s\in\{1\,,\dots,8\}$ and their corresponding estimators by $[\widehat{\boldsymbol{\beta}}]_s$\,.
For the one-layer model, we have
\begin{subequations}
	\label{eq:beta_1_4}
	\begin{gather}
		\underbrace{[\boldsymbol{\beta}]_1 = \left[\mkern-5mu\begin{array}{c}1500\\0.75\\0.0015\end{array}\mkern-5mu\right]}_{\text{near isotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_2 = \left[\mkern-5mu\begin{array}{c}1500\\0.75\\0.0408\end{array}\mkern-5mu\right]}_{\text{low anisotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_3 = \left[\mkern-5mu\begin{array}{c}1500\\0.75\\0.0832\end{array}\mkern-5mu\right]}_{\text{moderate anisotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_4 = \left[\mkern-5mu\begin{array}{c}1500\\0.75\\0.1728\end{array}\mkern-5mu\right]}_{\text{near isotropy}}
		\,;
		\tag{\theequation a,b,c,d}
	\end{gather}
\end{subequations}
for the two-layer model, we have
\begin{subequations}
	\label{eq:beta_5_8}
	\begin{gather}
		\underbrace{[\boldsymbol{\beta}]_5 = \left[\mkern-5mu\begin{array}{c} 911\\3285\\1.5\\0.5\\0.0015\\0.0019\end{array}\mkern-5mu\right]}_{\text{near isotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_6 = \left[\mkern-5mu\begin{array}{c} 911\\3285\\1.5\\0.5\\0.0408\\0.0618\end{array}\mkern-5mu\right]}_{\text{low anisotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_7 = \left[\mkern-5mu\begin{array}{c} 911\\3285\\1.5\\0.5\\0.0832\\0.1272\end{array}\mkern-5mu\right]}_{\text{moderate anisotropy}}
		\,,\quad
		\underbrace{[\boldsymbol{\beta}]_8 = \left[\mkern-5mu\begin{array}{c} 911\\3285\\1.5\\0.5\\0.1728\\0.2688\end{array}\mkern-5mu\right]}_{\text{high anisotropy}}
		\,.
		\tag{\theequation a,b,c,d}
	\end{gather}
\end{subequations}
As we recall from equation~\eqref{eq:chi}, the value of anisotropy corresponds to the difference between the horizontal and vertical velocity.
To quantify the amount of anisotropy, we set four categories based on the relative percentage increase in horizontal to vertical velocity.
They are: near isotropy, ranging from 0\% up to 1\%\,; low anisotropy, from 1\% up to 7\%\,; moderate anisotropy, from 7\% up to 15\%\,; and high anisotropy, from 15\% onward.
For model parameter sets~\eqref{eq:beta_1_4} and~\eqref{eq:beta_5_8}, the anisotropy category to which each belongs is indicated by its respective underbrace.

We restrict model-parameter sets to such values to remain consistent with reasonable values of anisotropy exhibited in shale~\citep[e.g.,][pp.~531--532]{VernikLiu1997}.
Finally, we tabulate the observed $[\widehat{\boldsymbol{\beta}}]_s$ and their corresponding $\delta([\widehat{\boldsymbol{\beta}}]_s)$ as well as discuss the reliability of results.

\subsection{One-layer optimization}
\label{sec:OneLayOpt}
For the one-layer model, we seek to solve problem~\eqref{eq:P(beta)}, for which
\begin{subequations}
	\label{eq:b>0,chi>0}
	\begin{gather}
		l_1(b_1) := \frac{1}{1+\exp\left(-10^6\,b_1\right)}
		\quad{\rm and}\quad
		l_2(\chi_1) := \frac{1}{1+\exp\left(-10^6\,\chi_1\right)}
		\,.
		\tag{\theequation a,b}
	\end{gather}
\end{subequations}
To do so, first, we ensure the convergence of the control experiment on model parameters~\eqref{eq:beta_1_4}, for which we require the use of Hessian modifications and logarithmic barrier functions, which are explained in Sections~\ref{sec:OptTech} and~\ref{sec:LogBar}, respectively, and justified in Appendix~\ref{app:N=1}.
Second, we repeat the optimization on model parameters~\eqref{eq:beta_1_4} under the influence of the six noise profiles plotted in Figure~\ref{fig:NoiseProfiles}.
We tabulate the resulting estimating vectors and their relative errors in Table~\ref{tab:OneLayOpt}.
Therein, we observe two prominent trends concerning the relative errors.

The first trend is the increase in relative error corresponding to the increase in the range of the noise profile.
We observe that for each increase in the order of magnitude of the noise profile, there is an increase by approximately one order of magnitude of the relative errors of the model parameters.
The second trend is the decrease in relative error corresponding to the increase in anisotropy.
We observe that, for the nearly isotropic case, the relative-error ordering is $\delta(a_1) < \delta(b_1) < \delta(\chi_1)$\,.
In contrast, as the amount of anisotropy increases, the relative error in $\chi_1$ diminishes to approximately the same level in $a_1$\,, with the largest relative error manifesting in $b_1$\,.
Thus, we state that, for the nearly isotropic case, the relative error manifests disproportionately in the anisotropy parameter whereas, for increasing levels of anisotropy, the relative error is spread out across all model parameters.
We expect this result because, for low anisotropy, small errors in $\chi_1$ greatly affect the modelled traveltimes, which is not the case for high anisotropy.
In other words, small errors in $\chi_1$ greatly affect the modelled traveltimes in low-anisotropy cases, but have less of an effect for increasing levels of anisotropy.

Along with these trends, we determine that the reliability of results is satisfactory for fixed-percentage noise ranges of at least~$\pm0.1\%$\,.
This percentage corresponds to noise-profile values on the order of milliseconds, which is the same order of magnitude associated with errors in picking traveltimes in real-data cases.
Thus, in the context of optimization, applying a one-layer model to traveltime data in regions where the subsurface is comprised predominantly of shale will yield reliable results.

\begin{table}[H]
	\centering
	\footnotesize
	\begin{tabular}{c*{6}{c}}
		\bottomrule
		& $\pm0.000001\%$ & $\pm0.00001\%$ & $\pm0.0001\%$ & $\pm0.001\%$ & $\pm0.01\%$ & $\pm0.1\%$ \\
		\toprule
		$[\widehat{\boldsymbol{\beta}}]_1$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.00005\\0.74999993\\0.00150001\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99990\\0.75000015\\0.00149998\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99866 \\ 0.75000198 \\ 0.00149968\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.98059 \\ 0.75002729 \\ 0.00149619\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.12664 \\ 0.74983824 \\ 0.00150612\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1498.25090 \\ 0.75254910 \\ 0.00109744\end{array}\mkern-10mu\right]$
		\\[12.5pt]
		$[\widehat{\boldsymbol{\beta}}]_2$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.00005 \\ 0.74999993 \\ 0.04080001\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99988 \\ 0.75000017 \\ 0.04079998\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99842 \\ 0.75000231 \\ 0.04079962\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.97821 \\ 0.75003055 \\ 0.04079562\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.13668 \\ 0.74982463 \\ 0.04080671\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1497.96923 \\ 0.75294017 \\ 0.04032501\end{array}\mkern-10mu\right]$
		\\[12.5pt]
		$[\widehat{\boldsymbol{\beta}}]_3$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.00006 \\ 0.74999992 \\ 0.08320001\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99987 \\ 0.75000020 \\ 0.08319998\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99815 \\ 0.75000269 \\ 0.08319955\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.97551 \\ 0.75003424 \\ 0.08319496\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.14766 \\ 0.74980972 \\ 0.08320740\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1497.64443 \\ 0.75339073 \\ 0.08263901\end{array}\mkern-10mu\right]$
		\\[12.5pt]
		$[\widehat{\boldsymbol{\beta}}]_4$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.00008 \\ 0.74999989 \\ 0.17280002\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99983 \\ 0.75000025 \\ 0.17279997\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.99749 \\ 0.75000360 \\ 0.17279937\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1499.96937 \\ 0.75004267 \\ 0.17279337\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1500.17137 \\ 0.74977744 \\ 0.17280899\end{array}\mkern-10mu\right]$ &
		$\left[\mkern-10mu\begin{array}{c} 1496.88716 \\  0.75443983 \\ 0.17202903\end{array}\mkern-10mu\right]$ \\
		\bottomrule
		& $\pm0.000001\%$ & $\pm0.00001\%$ & $\pm0.0001\%$ & $\pm0.001\%$ & $\pm0.01\%$ & $\pm0.1\%$ \\
		\toprule
		$\delta([\widehat{\boldsymbol{\beta}}]_1)$ &
		$\left[\mkern-10mu\begin{array}{r} 0.000003\\-0.000009\\0.000589\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000007\\0.000020\\-0.001186\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000089 \\ 0.000264 \\ -0.021257\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.001294\\ 0.003639 \\ -0.254192\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} 0.008443 \\ -0.021569 \\ 0.407760\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.116607 \\ 0.339880 \\ -26.83716\end{array}\mkern-5mu\right]$
		\\[12.5pt]
		$\delta([\widehat{\boldsymbol{\beta}}]_2)$ &
		$\left[\mkern-10mu\begin{array}{r}   0.000004 \\ -0.000010 \\ 0.000025\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000008 \\ 0.000023 \\ -0.000051\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.000105 \\ 0.000308 \\ -0.000931\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.001453 \\ 0.004073 \\ -0.010727\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}   0.009112 \\ -0.023382 \\ 0.016453\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.135382 \\ 0.392023 \\ -1.164195\end{array}\mkern-5mu\right]$
		\\[12.5pt]
		$\delta([\widehat{\boldsymbol{\beta}}]_3)$ &
		$\left[\mkern-10mu\begin{array}{r}  0.000004 \\ -0.000011 \\ 0.000015\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.000009 \\ 0.000026 \\ -0.000030\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000124 \\ 0.000359 \\ -0.000544\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.001633 \\ 0.004566 \\ -0.006055\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  0.009844 \\ -0.025371 \\ 0.008893\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.157038 \\ 0.452098 \\ -0.674271\end{array}\mkern-5mu\right]$
		\\[12.5pt]
		$\delta([\widehat{\boldsymbol{\beta}}]_4)$ &
		$\left[\mkern-10mu\begin{array}{r}  0.000005 \\ -0.000014 \\ 0.000001\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000011 \\ 0.000033 \\ -0.000019\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.000167 \\ 0.000480 \\ -0.000366\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} -0.002042 \\ 0.005690 \\ -0.003834\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r} 0.011425 \\ -0.029675 \\ 0.005205\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-10mu\begin{array}{r}  -0.207522 \\ 0.591977 \\ -0.446164\end{array}\mkern-5mu\right]$
		\\
		\bottomrule
	\end{tabular}
	\captionsetup{singlelinecheck=off}
	\caption[]{\small%
		Tabulated results of performing the inequality-constrained Newton's method in a one-layer model for $[\boldsymbol{\beta}]_1\,,\dots,[\boldsymbol{\beta}]_4$\,, for increasing ranges of constant-percentage noise profiles.
		}
	\label{tab:OneLayOpt}
\end{table}

\subsection{Two-layer optimization}
\label{sec:TwoLayOpt}
For the two-layer model, we seek to solve problem~\eqref{eq:P(beta)}, for which we require $b_i>0$ and $\chi_i>0$\,, for $i=1,2$\,.
As such, we require logarithmic barrier functions~\eqref{eq:b>0,chi>0}, along with another two for the second layer.
Then, we repeat the control experiment on model parameters~\eqref{eq:beta_5_8}, in the same manner described in Appendix~\ref{app:N=1}.
As with the one-layer case, we repeat the optimization under the influence of the six noise profiles plotted in Figure~\ref{fig:NoiseProfiles} and tabulate the resulting estimating vectors, and their relative errors, in Table~\ref{tab:TwoLayOpt}.
Therein, we observe the same trends concerning the relative errors described in Section~\ref{sec:OneLayOpt}---increasing noise-profile values result in increasing relative errors of model parameters, and increasing anisotropy results in greater error for all model parameters as opposed to anisotropy parameters only.

Pertaining to the relative errors in Table~\ref{tab:TwoLayOpt}, the threshold beyond which estimator vectors $\widehat{\boldsymbol{\beta}}$ cease to approximate the model parameters is between fixed-percentage ranges $\pm0.001\%$ and $\pm0.01\%$\,.
The corresponding modes of the orders of magnitude are microseconds and tens of microseconds, respectively, which are several orders of magnitude less than the one-layer-optimization threshold.
In contrast to the one-layer model, the increase in the number of model parameters from three to six reduces drastically the reliability of results for increasing ranges of constant-percentage noise profiles.
With reference to optimization, this is a substantial limitation in applying a two-layer model to real traveltime data.

Let us go further and consider a three-layer model; this would increase the number of model parameters to nine.
The increase from a one- to two-layer model reduces the reliability-of-results threshold from greater than the order of milliseconds to the order of microseconds.
We postulate that there would be---at least---a similar threshold reduction upon increasing the model complexity to three layers.
Apart from the fact that nine-parameter optimization may lead to many local minima and nonphysical results, the severe noise restrictions that beset the three-layer model render its optimization nonviable, especially under the influence of the amount of noise that may be present in real datasets, which we exemplify in Section~\ref{sec:RealDataExample}.
Thus, for pragmatic considerations, we do not consider the application of three-layer models to real traveltime data.

\begin{table}[H]
	\centering
	\scriptsize
	\begin{tabular}{c*{6}{c}}
		\bottomrule
		& $\pm0.000001\%$ & $\pm0.00001\%$ & $\pm0.0001\%$ & $\pm0.001\%$ & $\pm0.01\%$ & $\pm0.1\%$ \\
		\toprule
		$[\widehat{\boldsymbol{\beta}}]_5$ &
		$\left[\mkern-5mu\begin{array}{c} 911.009592\\3284.96820\\1.49998427\\0.50000299\\0.00149613\\0.00190835\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{c} 910.944258\\3285.07556\\1.50011945\\0.50001430\\0.00150593\\0.00187726\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{c} 910.266256\\3285.44169\\1.50170736\\0.50042727\\0.00150304\\0.00172162\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{c} 914.038158\\3273.10431\\1.49545520\\0.50203289\\0.00003179\\0.00493085\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{c} 905.603808\\3288.80029\\1.51216872\\0.50609521\\0.00205428\\0.00003065\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{c} 1036.74251\\3267.02054\\1.19449282\\0.56582609\\0.02503647\\0.00002969\end{array}\mkern-5mu\right]$
		\\
		$[\widehat{\boldsymbol{\beta}}]_6$ &
		$\left[\mkern-5mu\begin{array}{r} 911.012563 \\ 3284.96100 \\ 1.49997859 \\ 0.50000469 \\ 0.04079494 \\ 0.06181132\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 910.925996 \\ 3285.10399 \\ 1.50015794 \\ 0.50001434 \\ 0.04080930 \\ 0.06176587\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 910.039068 \\ 3285.71282 \\ 1.50220374 \\ 0.50049698 \\ 0.04082813 \\ 0.06150877\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 944.406832 \\ 3196.98077 \\ 1.44009001 \\ 0.51342586 \\ 0.02998654 \\ 0.08788399\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 839.871930\\3514.61164\\1.62708786\\0.47557539\\0.06823943\\0.00003125\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 977.211703\\3494.49641\\1.2799358\\0.53941197\\0.09791665\\0.00003108\end{array}\mkern-5mu\right]$
		\\
		$[\widehat{\boldsymbol{\beta}}]_7$ &
		$\left[\mkern-5mu\begin{array}{r} 911.016309 \\ 3284.95194 \\ 1.49997138 \\ 0.50000715 \\ 0.08319335 \\ 0.12721535\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 910.902569 \\ 3285.14096 \\ 1.50020731 \\ 0.50001289 \\ 0.08321418 \\ 0.12714955\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 909.745255 \\ 3286.08401 \\ 1.50284195 \\ 0.50056564 \\ 0.08326923 \\ 0.12673931\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 951.816141 \\ 3184.03462 \\ 1.42511862 \\ 0.51800472 \\ 0.06992063 \\ 0.16032199\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 786.378533 \\ 3753.38368 \\ 1.72096002 \\ 0.43061938 \\ 0.14168995 \\ 0.00003173\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 1032.35162 \\ 3461.13230 \\ 1.15858605 \\ 0.54271667 \\ 0.14797361 \\ 0.06828564\end{array}\mkern-5mu\right]$
		\\
		$[\widehat{\boldsymbol{\beta}}]_8$ &
		$\left[\mkern-5mu\begin{array}{r} 911.026614\\3284.92719\\1.49995141\\0.50001523\\0.17278860\\0.26882767\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 910.832797\\3285.25726\\1.50035339\\0.49999941\\0.17283187\\0.26869219\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 908.977690\\3287.06432\\1.50451227\\0.50068564\\0.17300230\\0.26780959\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 971.071974\\3151.18290\\1.38640710\\0.53100965\\0.15273713\\0.32179400\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 708.435786\\4260.41204\\1.85699375\\0.26060906\\0.30376789\\0.00003241\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} 1631.52423\\3212.87278\\0.05115720\\0.58844322\\0.23671151\\0.28866875\end{array}\mkern-5mu\right]$ \\
		\bottomrule
		& $\pm0.000001\%$ & $\pm0.00001\%$ & $\pm0.0001\%$ & $\pm0.001\%$ & $\pm0.01\%$ & $\pm0.1\%$ \\
		\toprule
		$\delta([\widehat{\boldsymbol{\beta}}]_5)$ &
		$\left[\mkern-2.5mu\begin{array}{r} 0.001053\\-0.000968\\-0.001048\\0.000598\\-0.257738\\0.439608\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.006119\\0.002300\\0.007963\\0.002859\\0.395071\\-1.197096\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.080543\\0.013446\\0.113824\\0.085455\\0.202615\\-9.388697\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 0.333497\\-0.362122\\-0.302987\\0.406579\\-97.88077\\159.5187\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.592337\\0.115686\\0.811248\\1.219041\\36.95199\\-98.38681\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 13.80269\\-0.547320\\-20.36715\\13.16522\\1569.098\\-98.43740\end{array}\mkern-2.5mu\right]$
		\\
		$\delta([\widehat{\boldsymbol{\beta}}]_6)$ &
		$\left[\mkern-2.5mu\begin{array}{r} 0.001379 \\ -0.001187 \\ -0.001427 \\ 0.000937 \\ -0.012400 \\ 0.018309\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.008123 \\ 0.003166 \\ 0.010529 \\ 0.002867 \\ 0.022791 \\ -0.055229\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-5mu\begin{array}{r} -0.105481 \\ 0.021699 \\ 0.146916 \\ 0.099395 \\ 0.068935 \\ -0.471239\end{array}\mkern-5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 3.667051 \\ -2.679429 \\ -3.993999 \\ 2.685172 \\ -26.50357 \\ 42.20711\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -7.807692\\6.989700\\8.472524\\-4.884922\\67.25350\\-99.94944\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 7.268025\\6.377364\\-14.67095\\7.882394\\139.9918\\-99.94970\end{array}\mkern-2.5mu\right]$
		\\
		$\delta([\widehat{\boldsymbol{\beta}}]_7)$ &
		$\left[\mkern-2.5mu\begin{array}{r} 0.001790 \\ -0.001463 \\ -0.001908 \\ 0.001429 \\ -0.007997 \\ 0.012063\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.010695 \\ 0.004291 \\ 0.013820 \\ 0.002578 \\ 0.017049 \\ -0.039664\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.137733\\0.032999\\0.189463\\0.113128\\0.083210\\-0.362174\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 4.480367 \\ -3.073527 \\ -4.992092 \\ 3.600944 \\ -15.96078 \\ 26.03930\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -13.67963 \\ 14.25826 \\ 14.73067 \\ -13.87612 \\ 70.30042 \\ -99.97505\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 13.32071 \\ 5.361714 \\ -22.76093 \\ 8.543335 \\ 77.85290 \\ -46.31632\end{array}\mkern-2.5mu\right]$
		\\
		$\delta([\widehat{\boldsymbol{\beta}}]_8)$ &
		$\left[\mkern-2.5mu\begin{array}{r} 0.002921\\-0.002216\\-0.003239\\0.003047\\-0.006598\\0.010294\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.018354\\0.007831\\0.023560\\-0.000118\\0.018446\\-0.040109\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -0.221988\\0.062841\\0.300818\\0.137128\\0.117069\\-0.368455\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 6.594070\\-4.073580\\-7.572860\\6.201931\\-11.61045\\19.71503\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} -22.23537\\29.69291\\23.79958\\-47.87819\\75.79161\\-99.98794\end{array}\mkern-2.5mu\right]$ &
		$\left[\mkern-2.5mu\begin{array}{r} 79.09157\\-2.195654\\-96.58952\\17.68864\\36.98583\\7.391650\end{array}\mkern-2.5mu\right]$ \\
		\bottomrule
	\end{tabular}
	\captionsetup{singlelinecheck=off}
	\caption[]{\small%
		Tabulated results of performing the inequality-constrained Newton's method in a two-layer model for $[\boldsymbol{\beta}]_5\,,\dots,[\boldsymbol{\beta}]_8$\,, for increasing ranges of constant-percentage noise profiles.
	}
	\label{tab:TwoLayOpt}
\end{table}

\section{Real-data example}
\label{sec:RealDataExample}
Let us perform the optimization on a set of measured traveltimes, $\boldsymbol{T}$\,, obtained in offshore Newfoundland~\citep{Kaderali2009}.
The purpose is to use traveltime model~\eqref{eq:t} to solve problem~\eqref{eq:P(beta)} and, in turn, quantify the anisotropy and inhomogeneity of the subsurface in the region.

The subsurface unit in question is comprised predominantly of shale units that can be amalgamated effectively to a single unit.
We are able to corroborate this statement by considering the interval velocities corresponding to the dataset, which we plot in Figure~\ref{fig:Checkshot}.
Therein, we observe a steady linear increase in velocity throughout the entirety of the 2\,km subsurface.
As such, for this dataset, we consider a one-layer model to be an adequate, albeit low-resolution, representation of the subsurface.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.25\textwidth]{Fig_Checkshot.pdf}
	\caption{\small Interval velocities with depth}
	\label{fig:Checkshot}
\end{figure}

Further inspection of Figure~\ref{fig:Checkshot} reveals that there are zones of changing interval velocities, namely at depths $z=1212$\,m and $z=1585$\,m\,.
Upon this basis, it is reasonable to apply a multilayer model that segments the subsurface, depending on the desired level of resolution.
Thus, for the two-layer model, we bisect the subsurface and set the layer interface at depth $z_1 = 1212$\,m\,.
We select this layer-interface depth instead of $z=1585$\,m to avoid creating a model with a single layer that dominates the effect of the modelled traveltimes.

Prior to performing the optimizations, we must assess the amount of noise present in the real data.
To this end, we apply a three-point median filter on the traveltimes and calculate the difference between $\boldsymbol{T}$ and the filtered traveltimes to obtain the noise profile.
We plot the noise profile in Figure~\ref{fig:tReal_3pt_NoiseProfile} and, through visual inspection, we conclude that the amount of noise is on the order of milliseconds.
We calculate that the mode of the orders of magnitude of the noise profile is $10^{-4}$\,, which indicates that the most frequent noise-profile value is on the order of a fraction of a millisecond.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.5\textwidth]{Fig_tReal_3pt_NoiseProfile.pdf}
	\caption{\small Noise profile of measured traveltimes $\boldsymbol{T}$ obtained using a three-point median filter}
	\label{fig:tReal_3pt_NoiseProfile}
\end{figure}

As we recall from Sections~\ref{sec:OneLayOpt} and~\ref{sec:TwoLayOpt}, the reliability-of-results threshold for one-layer optimization is greater than milliseconds and for two-layer optimization is on the order of microseconds.
As such, we recognize that, in the context of the optimization, applying a one-layer model to traveltimes $\boldsymbol{T}$ will yield reliable results, but not so for the two-layer model.
Be that as it may, we perform the optimization for both one- and two-layer models, with restrictions $0<b_i<2$ and $0<\chi_i<1$\,.
The estimator vector and objective-function value, for the former, are
\begin{subequations}
	\label{eq:hat{beta}_N=1,hat{f}_N=1}
	\begin{gather}
		[\widehat{\boldsymbol{\beta}}]_{N=1} = \left[\mkern-2.5mu\begin{array}{r} 1293.79163 \\ 1.12823361 \\ 0.05416071 \end{array}\mkern-2.5mu\right]
		\quad{\rm and}\quad
		\hat{f}_{N=1} = 0.00023430
		\,,
		\tag{\theequation a,b}
	\end{gather}
\end{subequations}
and, for the latter, are
\begin{subequations}
	\label{eq:hat{beta}_N=2,hat{f}_N=2}
	\begin{gather}
		[\widehat{\boldsymbol{\beta}}]_{N=2} = \left[\mkern-2.5mu\begin{array}{r} 897.950408 \\ 3296.42981 \\ 1.98094661 \\ 0.00003242 \\ 0.00002178 \\ 0.00002216\end{array}\mkern-2.5mu\right]
		\quad{\rm and}\quad
		\hat{f}_{N=2} = 0.00009644
		\,.
		\tag{\theequation a,b}
	\end{gather}
\end{subequations}
To gain insight on the values of objective functions~(\ref{eq:hat{beta}_N=1,hat{f}_N=1}b) and~(\ref{eq:hat{beta}_N=2,hat{f}_N=2}b), let us use model parameters~(\ref{eq:hat{beta}_N=1,hat{f}_N=1}a) and~(\ref{eq:hat{beta}_N=2,hat{f}_N=2}a) to calculate the best-fit traveltimes using equation~\eqref{eq:t}.
Taking their difference with traveltimes $\boldsymbol{T}$ for every offset, we obtain the traveltime residuals for  both optimizations, which we plot in Figure~\ref{fig:tReal_Res}.
Therein, we observe that both traveltime residuals vary with offset.
Since the general trends of both the residuals remain within a $\pm$2\,ms range, we maintain that both fits are sufficiently within the traveltime picking error.
However, it is clear that both residual plots exhibit a cyclic pattern, with Figure~\ref{fig:tReal_Res_1Lay}, the one-layer case, being more pronounced than Figure~\ref{fig:tReal_Res_2Lay}, the two-layer case.
We attribute the cyclic pattern to an inadequacy of the model to simulate the measured traveltimes.
The inadequacy may be rooted in an insufficient number of layers to model the subsurface or an absence of elliptical velocity dependence in all layers.
Within the context of optimization, we expect the residuals to be random and appear to be centred about zero, which is the case of the noise profile of the measured traveltimes in Figure~\ref{fig:tReal_3pt_NoiseProfile}.
To this end, increasing the number of layers not only improves the representation of the medium, but also reduces the cyclic pattern, as can be observed in Figure~\ref{fig:tReal_Res_2Lay}.
With that being said, the cyclic pattern persists, which could be the consequence of imposing elliptical anisotropy in all layers.
Let us examine Figures~\ref{fig:tReal_Res_1Lay} and~\ref{fig:tReal_Res_2Lay} individually.

\captionsetup[subfigure]{subrefformat=simple,labelformat=simple}
\renewcommand\thesubfigure{(\alph{subfigure})}
\begin{figure}[h]
	\centering
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Fig_tReal_Res_1Lay.pdf}
		\caption{}
		\label{fig:tReal_Res_1Lay}
	\end{subfigure}%
	~
	\begin{subfigure}[t]{0.49\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Fig_tReal_Res_2Lay.pdf}
		\caption{}
		\label{fig:tReal_Res_2Lay}
	\end{subfigure}
	\caption{\small%
		Traveltime residuals with offset upon performing the (a) one- and (b) two-layer optimization on traveltimes~$\boldsymbol{T}$
	}
	\label{fig:tReal_Res}
\end{figure}

The traveltime residual in Figure~\ref{fig:tReal_Res_1Lay} does not remain near zero with any offset, which could be attributed to a suboptimal model choice.
Indeed, consolidating nearly two kilometres of subsurface into a single layer does not necessarily lead to an accurate representation of its anisotropy and inhomogeneity.
However, since the subsurface is considered to be predominantly shale and the general trend of the traveltime residual remains within the experimental picking error limits, we claim that this model is an adequate representation of the subsurface, at least at a macroscopic level.
Furthermore, in view of the one-layer-optimization threshold described in Section~\ref{sec:OneLayOpt}, the amount of noise exhibited in $\boldsymbol{T}$ should not lead to poor numerical results.
Thus, we contend that model parameters~(\ref{eq:hat{beta}_N=1,hat{f}_N=1}a) are reliable.

Recalling the adage that an increase in the number of model parameters results in an increase of numerical fit, the traveltime residual in Figure~\ref{fig:tReal_Res_2Lay} remains closer to zero than in Figure~\ref{fig:tReal_Res_1Lay}, and hence, $\hat{f}_{N=2} < \hat{f}_{N=1}$\,, as expected.
Indeed, bisecting the subsurface improves the numerical fit of traveltimes $\boldsymbol{T}$\,.
However, in view of the reliability-of-results threshold for two-layer optimization described in Section~\ref{sec:TwoLayOpt}, the amount of noise exhibited in $\boldsymbol{T}$\,, which is on the order of milliseconds, is three orders of magnitude greater than the threshold, which is on the order of microseconds.
As such, the reliability of model parameters~(\ref{eq:hat{beta}_N=1,hat{f}_N=1}b) is poor; they indicate that the first layer is strongly inhomogeneous and exhibits low anisotropy, and that the second layer is nearly homogeneous and exhibits near isotropy.
Such results are not consistent with the interval velocities plotted in Figure~\ref{fig:Checkshot}, nor with the knowledge that the subsurface in question is comprised predominantly of shale.
Thus, despite the increase of model adequacy for physical representation of the subsurface and the numerical fit of traveltimes, there is a significant tradeoff in the reliability of results.

\section{Conclusions}
To estimate the values of parameters, attributed to anisotropy and inhomogeneity within the context of an analytical mathematical model, using vertical seismic profile traveltimes, we use Newton's method to minimize the least-squares residual between measured and modelled traveltimes.
To understand applications of this method to real traveltimes, we construct synthetic models to analyze the reliability of results under the influence of increasing percentage ranges of noise.
We determine that there is a threshold for the acceptable level of noise, which increases in severity as the model complexity increases.
In particular, the reliability-of-results threshold for a one-layer model is three orders of magnitude greater than that of a two-layer model.

Within the context of traveltime equation~\eqref{eq:t}, greater complexity is tantamount to segmenting the subsurface into a greater number of layers.
Insofar as constructing physically representative models of the subsurface is concerned, the greater the model complexity, the better the fit, as expected due to the larger number of parameters used to accommodate the data.
However, as determined in Section~\ref{sec:RealDataExample}, the amount of noise present in real traveltime data imposes a tradeoff in the reliability of results with an increase in model complexity.

Pertaining to the real-data example considered in Section~\ref{sec:RealDataExample}, we find that, despite our desire to construct a multilayer model, the noise present in this dataset exceeds the reliability-of-results threshold  for a two-layer model.
Although we are able to minimize the traveltime residual for a model comprised of two or more layers, we are restricted to a one-layer model for reliable results.
However, since the subsurface in question is comprised predominantly of shale, we contend that the one-layer model is an adequate representation and that model parameters~\eqref{eq:hat{beta}_N=1,hat{f}_N=1} are reliable and acceptable.
\section*{Acknowledgments}
We wish to acknowledge discussions with Michael A. Slawinski and the graphic support of Elena Patarini.
Also, this research was partially supported by the Natural Sciences and Engineering Research Council of Canada, grant 202259.
\bibliographystyle{apa}


\section{Introduction}
\label{sec1}
Day-ahead electricity markets are multi-unit auctions where market players submit orders to buy or sell electric power in a given area during predefined hours of the following day. Market players agree on a collection of rules describing the bid acceptance and the price determination. The so called market operator collects the participants' orders, and determines the dispatch in order to maximize the total social welfare (TSW).

The TSW is defined as the total utility of consumption minus the total cost of production, and it may be derived as a function of the decision variables describing the acceptance/rejection of bids and bid parameters. Assuming a general market clearing price, the total social welfare may be decomposed to the welfare or surplus of single bids.

The computed market prices should ideally support a market equilibrium, which maximizes the social welfare \cite{madani2017revisiting}, with respect to the actual constraints (e.g. keeping production and consumption in balance for every hour). \\
The most difficult feature to deal with in day-ahead electricity markets is the fact, that some orders may be non-convex \cite{madani2016non}.

Non-convexity may originate from technical constraints such as minimum performance or load gradient \cite{sleisz2019new}, and can be also generated by the cost structure (e.g. start-up costs of power plants). \\
Regarding non-convex bids, in this study we assume the presence of block orders (orders including potentially multiple bids for several consecutive hours \cite{meeus2009block}) for which a "fill-or-kill condition" must met: the order is either fully accepted or fully rejected (in each of the hours).  These block orders help participants to express their production constraints and cost structures: they allow the participant to combine multiple periods.

\section{Methods}
The input of the day-ahead market is the set of supply and demand bids for each period (hour). Each submitted bid is characterized by a quantity and a bid price (per unit). It is our task to decide which bids will be accepted and rejected.  Bids which define no interdependencies between periods and may be partially accepted are called standard (hourly) bids.

\subsection{General multiperiod problem formulation with block bids}
The general problem considered in this paper is including $T$ time periods (hours), (multi-period) block bids, minimum allocation ratios and start-up costs as well.
Each block bid $j$ is characterized by a set $T_j \subseteq T$, describing the (consecutive) time periods for which the block bid is submitted, by a quantity $Q_j$ and by a bid price $P_j$. We assume that the latter two parameters are universal for all periods included in the bid (in other words we do not consider 'profiled block bids' \cite{euphemia2015}).

\begin{align}
& \max\limits_{x, u} \sum_{t \in T} \left( \sum_{i \in I_t} Q_{ti}P_{ti} x_{ti} + \sum_{j: t\in T_j} Q_j P_j u_j \right) - \nonumber \\ & - \sum_{t \in T}\sum_{i \in I_t}u_{ti}F_{ti} \label{obj fncn} \\
\text{wrt.} ~~~~ &r_{ti}u_{ti} \leq x_{ti}\leq u_{ti}~~~~~~ \forall t \in T, ~~~\forall i \in I_t \label{startup}\\
 &Q_{ti} x_{ti} + \sum_{j: t\in T_j} Q_j u_j = 0~~~~\forall t \in T \label{power balance}
\end{align}
where $x_{ti} \in [0,1]$ denotes the acceptance indicator of the $i$-th standard hourly bid submitted for period $t$, $r_{ti}$ denotes its minimum allocation ratio and $I_t$ is the index set of standard hourly bids submitted for period $t$.
$u_j \in \{0,1 \}$ denotes the binary acceptance indicator of block bid $j$.
The first term in the objective function (\ref{obj fncn}) corresponds to the TSW of hourly bids while the second term corresponds to the welfare contribution of block bids. $F_{ti}$ corresponds to the start-up cost for the $i$-th standard hourly bid of period $t$. If $F_{ti}=0$, the corresponding variable $u_{ti}$ may be fixed to 1, and it may be neglected in the optimization problem.
We assume that in the case of block bids, bidders include the start-up costs in the bid price as well.
Equation (\ref{power balance}) formulates the power balance for all hours, while eq. (\ref{startup})
corresponds to the consideration that without the start-up of the respective unit, no bid may be accepted.

\subsection{The concept of bid-aggregation based clearing}
The number of bids in the clearing problem naturally affects the computational time of the optimization.
By reducing the number of bids, the speed of optimization can be increased.

We introduce the concept of bid aggregation to decrease the number of bids (thus accelerate the clearing process).
Solving the aggregated optimization problem, based on the  acceptance/rejection of certain aggregated bids we will fix the acceptance indicators of the corresponding original bids, and we may exclude them from the further computational steps aiming the solution of the original problem.
In these further steps we decompose the remaining (not fixed) aggregated bids to determine the acceptance of the remaining original bids. Multi-period block bids are not included in the aggregation process.

\subsubsection{The aggregation rule for $n$ bids}
\label{aggreg_rule_for_n_bids}
$n$ standard supply/demand offers are given that we want to aggregate into a single bid. The parameters\footnote{$P^A$ stands for the unit price, $Q^A$ for the quantity, $r^A$ for the minimum allocation ratio and $F^A$ for the start-up cost of the aggregated bid.} of the resulting aggregated bid will be computed as follows:
\begin{equation*}
\begin{aligned}
&P^A = \frac{\sum_{i = 1}^{n} p_i q_i}{\sum_{i = 1}^{n} q_i} ~~~~Q^A = \sum_{i = 1}^{n} q_i \nonumber \\
&r^A = \frac{\sum_{i = 1}^{n} r_i q_i}{\sum_{i = 1}^{n} q_i} ~~~~F^A = \sum_{i = 1}^{n} F_i \label{aggregation}
\end{aligned}
\end{equation*}

The bid price of the aggregated bid ($P^A$) is computed from the quantity-weighted average of original bid prices, while its quantity ($Q^A$) will be the sum of the quantities of its bids. The minimum allocation ratio of the aggregated bid ($r^A$) is the quantity-weighted average of the minimum allocation ratios of its components and its startup cost ($F^A$) is determined as the sum of the startup costs of the original bids.

\subsection{Clearing of the aggregated bids - the first step of the clearing}
\label{Cotab}
By partitioning the original bids into subsets and aggregating these subsets into new bids we decrease the size of the optimization problem.
The solution of this new optimization problem results in the determination of acceptance indicators for the aggregated bids.

In most of the cases
there are three possible outcomes of this first optimization step:
 \begin{itemize}
   \item One of the demand bids is partially accepted.
   \item One of the supply bids is partially accepted.
   \item The acceptance indicator of all bids is 0 or 1.
 \end{itemize}

If there is an aggregated bid on the supply (demand) side that has been partially accepted, we can identify the two closest bids on the demand (supply) side: The first fully accepted demand (supply) bid whose price is higher than the bid price of the partially accepted bid, and the first fully rejected demand (supply) bid  whose price is lower than the bid price of the partially accepted bid.

Thus, we have three distinguished aggregated bids: one partially accepted (on one side: supply or demand), one fully accepted and one fully rejected (the latter two on the other side).

If the acceptance indicator of all bids is 0 or 1 we take four aggregated bids similarly as described above (one-one fully accepted and rejected bids from both the supply and demand side).

\subsection{The second clearing step}
The initial bids from which these three (or four) distinguished aggregated bids are composed may be identified by the stored IDs. In the next step of the process we first fix the acceptance value of all original bids, which are not included in these three (or four) distinguished aggregated bids according to the outcome of the first clearing step (to fully accepted or rejected), and perform a second optimization step to determine the final acceptance values of those bids which are included in the three (or four) distinguished aggregated bids.

\subsubsection{Constraints of the second clearing step}
Let us remember that in the original problem eq. (\ref{power balance}) ensures the balance of supply and demand quantities for each period. This balance constraints also hold in the first step of the process (when clearing the aggregated bids).

However, if according to the outcome of this first step we fix the acceptance indicators of some bids (potentially both on the supply and on the demand side), in most cases we create a power imbalance. In order to get a final result which respects the power balance constraint, during the second step of the process (in the second clearing step), this imbalance must be taken into account and the remaining bids must be cleared in a way which compensates for this imbalance. This can be done via the appropriate modification of the right-hand side of eq. (\ref{power balance}).

\subsection{Example for the aggregation-based clearing}
\label{example}
Let us consider a small one-period example given in Tables \ref{ex_aggr_dem} and \ref{ex_aggr_sup} without minimum allocation ratios, startup costs and block bids for the demonstration of the bid-aggregation based clearing:

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
 Bid ID & $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$  \\ \hline
 1 & $23$ & $78$ \\ \hline
 2 & $67$ & $65$ \\ \hline
 3 & $27$ & $57$ \\ \hline
 4 & $30$ & $55$ \\ \hline
 5 & $91$ & $45$ \\ \hline
 6 & $90$ & $42$ \\ \hline
\end{tabular}
\caption{Demand offers \label{ex_aggr_dem}}%
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
Bid ID & $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$  \\ \hline
 7 & $-96$ & $40$ \\ \hline
 8 & $-71$ & $42$ \\ \hline
 9 & $-41$ & $47$ \\ \hline
 10 & $-80$ & $52$ \\ \hline
 11 & $-99$ & $57$ \\ \hline
 12 & $-99$ & $58$ \\ \hline
\end{tabular}
\caption{Supply offers \label{ex_aggr_sup}}%
\end{table}

The result of solving the problem without aggregation can be seen in Table \ref{res_wo_aggr}.
\begin{table}[h!]
\centering
\caption{Results without aggregation \label{res_wo_aggr}}
\begin{tabular}[t]{ |c|>{\centering\arraybackslash}p{1cm}||c|>{\centering\arraybackslash}p{1cm}|}
 \hline
Bid ID &$x_i$ & Bid ID &$x_i$ \\ \hline
 1 & $1$ & 7 & $1$ \\ \hline
 2 & $1$ & 8 & $1$ \\ \hline
 3 & $1$ & 9 & $0$ \\ \hline
 4 & $1$ & 10 &$0$ \\ \hline
 5 & $0.2198$ & 11 &$0$ \\ \hline
 6 & $0$ & 12 &$0$ \\ \hline
\end{tabular}
\end{table}

Let us now consider pairwise aggregation of the original bids. The IDs of the created aggregated bids are denoted by A$i$. In this case, according to the formulas defined in  (\ref{aggregation}) the aggregated bids are as follows (Tables \ref{aggr_dem} and \ref{aggr_sup}).

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
 Aggr. bid ID& $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$ & Init. bids \\ \hline
 A1 & $90$ & $68.32$ & 1, 2 \\ \hline
 A2 & $57$ & $55.95$ & 3, 4 \\ \hline
 A3 & $181$ & $43.51$ & 5, 6 \\ \hline
\end{tabular}
\caption{Aggregated demand offers \label{aggr_dem}}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
 Aggr. bid ID& $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$ & Init. bids \\ \hline
 A4 & $-167$ & $40.85$ & 7, 8 \\ \hline
 A5 & $-121$ & $50.31$ & 9, 10 \\ \hline
 A6 & $-198$ & $57.5$ & 11, 12 \\ \hline
\end{tabular}
\caption{Aggregated supply offers \label{aggr_sup}}
\end{table}

Figures \ref{fig:rendes_licitek} and \ref{fig:aggregalt_licitek} depict the supply-demand curves of the original problem and of the problem derived from the aggregated bids.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.4]{rendes_licitek.png}
\caption{Initial bids \label{fig:rendes_licitek}}%
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[scale=0.4]{aggregalt_licitek.png}
\caption{Aggregated bids \label{fig:aggregalt_licitek}}%
\end{figure}

Solving the optimization problem described in eqs. (\ref{obj fncn}) - (\ref{power balance}) for the aggregated bids we get:

\begin{table}[h!]
\centering
\caption{Resulting acceptance indicators of the aggregated bids ($x^A_i$) \label{aggr_ind}}
\begin{tabular}{ |c|>{\centering\arraybackslash}p{1cm}||c|>{\centering\arraybackslash}p{1cm}| }
 \hline
Aggr. bid ID &$x_i^A$ & Aggr. bid ID &$x_i^A$ \\ \hline
 A1 & $1$ & A4 & $1$ \\ \hline
 A2 & $1$ & A5 & $0$ \\ \hline
 A3 & $0.1105$ & A6 & $0$ \\ \hline
\end{tabular}
\end{table}

As can be seen in Figure \ref{fig:aggregalt_licitek} we partially accept the A3 aggregated  demand bid.
Thus, the set of the three distinguished aggregated bids defined in subsection \ref{Cotab} will be A3, A4 and A5 (see Figure \ref{fig:szetbontando_aggregalt_licitek}) --  the acceptance indicators of the components of these bids (regarding original bid IDs: 5,6,7,8,9,10) will serve as variables of the second clearing step, while the acceptance of all other original bids are fixed according to the resulting acceptance of the aggregated bid including them as summarized in Table \ref{fixed_ind}.

\begin{table}[h!]
\begin{center}
\begin{tabular}{ |c|c|c| }
 \hline
Bid ID &$x_i$ & based on\\ \hline
 1 & $1$ & full acceptance of $x^A_1$\\ \hline
 2 & $1$ & full acceptance of $x^A_1$\\ \hline
 3 & $1$ & full acceptance of $x^A_2$\\ \hline
 4 & $1$ & full acceptance of $x^A_2$\\ \hline
 11 & $0$ & full rejection of $x^A_6$\\ \hline
 12 & $0$ & full rejection of $x^A_6$\\ \hline
\end{tabular}
\caption{ Fixed variables regarding the acceptance indicators of the original problem. \label{fixed_ind}}
\end{center}
\end{table}

After this, we perform the second optimization with the remaining original bids (Tables \ref{remaining_dem} and \ref{remaining_sup}) can be seen in Figure \ref{fig:szetbontott_aggregalt_licitek}.

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
 Bid ID & $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$  \\ \hline
 5 & $91$ & $45$ \\ \hline
 6 & $90$ & $42$ \\ \hline
\end{tabular}
  \caption{Remaining demand offers \label{remaining_dem}}
\end{table}

\begin{table}[h!]
\centering
\begin{tabular}{ |c|c|c|c|c| }
 \hline
 Bid ID & $Q\> [MWh]$ & $P\> [\text{EUR} / MWh]$  \\ \hline
 7 & $-96$ & $40$ \\ \hline
 8 & $-71$ & $42$ \\ \hline
 9 & $-41$ & $47$ \\ \hline
 10 & $-80$ & $52$ \\ \hline
\end{tabular}
  \caption{Remaining supply offers \label{remaining_sup}}
\end{table}

As mentioned earlier, fixing the acceptance indicators of bids not corresponding to the three distinguished aggregated bid creates a power imbalance, which must be compensated in the second optimization step. In this example we fix the acceptance indicators of bids 1, 2, 3 and 4 to 1 (see Table \ref{aggr_ind}), resulting in a surplus of
\begin{equation*}
q_1 + q_2 + q_3 + q_4 = 23 + 67 + 27 + 30 = 147 MWh
\end{equation*}
in the demand side, and this has to be taken into account.

The value required to satisfy the power balance constraint calculated above is placed in the appropriate position of the $b_{eq}$ vector (its position in the vector is the number of the observed period).

As we need to compensate the demand surplus with production, and supply bids do have a negative sign for quantity, the balance equation for the second optimization is as follows:
\begin{equation*}
\begin{aligned}
&  91 \cdot x_5 + 90 \cdot x_6 -96 \cdot x_7 - 71 \cdot x_8 - \\ & - 41  \cdot x_9 -80 \cdot x_{10} = -147.
\end{aligned}
\end{equation*}

Figures \ref{fig:szetbontando_aggregalt_licitek} and \ref{fig:szetbontott_aggregalt_licitek} depict the remaining aggregated and decomposed offers.
\begin{figure}[h!]
    \centering
\includegraphics[scale=0.4]{szetbontando_aggregalt_licitek.png}
  \caption{The 3 distinguished aggregated bids, which will be decomposed in the second step.\label{fig:szetbontando_aggregalt_licitek}}
\end{figure}

\begin{figure}[h!]
    \centering
\includegraphics[scale=0.4]{szetbontott_aggregalt_licitek.png}
  \caption{The remaining original bids after the decomposition of the 3 distinguished aggregated bids. \label{fig:szetbontott_aggregalt_licitek}}
\end{figure}

Solving the second clearing step gives the result (see Figure \ref{fig:szetbontott_aggregalt_licitek}):

\begin{table}[h!]
\centering
\caption{Bids indicators' in the second clearing step \label{sec_opt_ind}}
\begin{tabular}{ |c|>{\centering\arraybackslash}p{1cm}||c|>{\centering\arraybackslash}p{1cm}| }
 \hline
Bid ID &$x_i$ & Bid ID &$x_i$\\ \hline
 5 & $0.2198$ &  8 & $1$\\ \hline
 6 & $0$ &  9 &$0$ \\ \hline
 7 & $1$ &  10 &$0$ \\ \hline
\end{tabular}
\end{table}

Testing the quantities we get that
\begin{equation*}
\begin{aligned}
& x_5 \cdot q_5 + x_6 \cdot q_6 - x_7 \cdot q_7 - \\ & - x_8 \cdot q_8 - x_9 \cdot q_9 - x_{10} \cdot q_{10} = \\
&= 0.2198 \cdot 91 + 0 \cdot 90 - 1 \cdot 96 - \\ & - 1 \cdot 71 - 0 \cdot 41 - 0 \cdot 80 = - 147,
\end{aligned}
\end{equation*}
the requirement for the power balance is satisfied.\\

The final result using the results from the first (see Table \ref{fixed_ind}) and the second (see Table \ref{sec_opt_ind}) steps of optimization:

\begin{table}[h!]
\centering
\caption{Results with aggregation\label{result}}
\begin{tabular}[t]{ |c|>{\centering\arraybackslash}p{1cm}||c|>{\centering\arraybackslash}p{1cm}|}
 \hline
Bid ID &$x_i$ & Bid ID &$x_i$ \\ \hline
 1 & $1$ & 7 & $1$ \\ \hline
 2 & $1$ & 8 & $1$ \\ \hline
 3 & $1$ & 9 & $0$ \\ \hline
 4 & $1$ & 10 &$0$ \\ \hline
 5 & $0.2198$ & 11 &$0$ \\ \hline
 6 & $0$ & 12 &$0$ \\ \hline
\end{tabular}
\end{table}

In this case, the results of the original method and the method based on the aggregation of bids are the same. In the following however, we will see that the bid aggregation-based clearing may also induce sub-optimal results in some cases.

\section{Results and discussion}
In the previous example (\ref{example}), we have always grouped 2 bids into an aggregated bid (thus $n$ was equal to 2), but for larger problems, aggregating a low number of bids does not results in the most efficient solution because the time needed for the first optimization step will be higher due to the large number of aggregated bids.

If we perform the aggregation for flexible number of bids, taking into account their price parameters, we potentially get an aggregated bid set with more equidistant distribution in the price space, which is beneficial for optimization purposes. We will call the proposed method \emph{nominal aggregation}.

Although in the majority of the tested scenarios the nominal aggregation based clearing reaches the same TSW as the non-aggregated clearing (thus the result is optimal), in some cases a sub-optimal solution is reached.
It is possible, that the decomposition of aggregated bids not fixed after the first clearing step results in a bid set, which can not lead to an optimal solution.

Figures \ref{fig:nominalis_nem_talalja_el_3} and \ref{fig:nominalis_nem_talalja_el_4} illustrate such a case: The
remaining 3 aggregated bids and the initial bids contained by them are depicted in the figures.
\begin{figure}[h!]
\centering
\includegraphics[scale=0.5]{nominalis_nem_talalja_el_3.eps}
\caption{The 3 distinguished aggregated bids, which will be decomposed in the second step.
\label{fig:nominalis_nem_talalja_el_3}}
\end{figure}

\begin{figure}[h!]
\centering
\includegraphics[scale=0.5]{nominalis_nem_talalja_el_4.eps}
\caption{The remaining original bids after the decomposition of the 3 distinguished aggregated bids. \label{fig:nominalis_nem_talalja_el_4}}
\end{figure}

As can be seen on Fig. \ref{fig:nominalis_nem_talalja_el_3} the curve of the aggregated supply bid meets with the aggregated demand curve (this is the reason why these three aggregated bids were chosen to participate in the second optimization).
However in Fig. \ref{fig:nominalis_nem_talalja_el_4} the supply and demand curves after the decomposition do not intersect each other.
If we do not use bid aggregation and solve the original problem, in the TSW maximizing solution, bid 13 is accepted partially.
In contrast, using bid aggregation, its acceptance indicator will be fixed to 0 (according to the acceptance indicators of aggregated offers). From this, it can be clearly seen that the optimal solution can not be reached here via the bid aggregation method.

\subsection{Optimality performance}
Figures \ref{results1} and \ref{results2} depict the suboptimality of the aggregation-based method with the number of bids on the $x$-axis. We examined the results in 250 randomly generated test cases using a predetermined threshold ($\varepsilon$), assuming various number of bids and varying penetration of block bids in the market. If the relative deviation of the TSW value obtained by the actual aggregation method was less then $\varepsilon$, we considered the clearing successful. The maximum of the nominal aggregation's average relative difference from the reference values was 0.3422\% (It occured when the number of supply-demand bids was 100 and the ratio of block bids was 20\%.).

\begin{figure}[h!]
    \centering
\begin{tikzpicture}
\pgfplotsset{height=7cm, width=\textwidth,compat=1.17}
\begin{axis}[
ybar,
enlarge y limits={0.2},
enlarge x limits={0.5},
legend style={at={(0.5,1), anchor=north west},
anchor=north,legend columns=-1},
ylabel={Hit rate [\%]},
xlabel={Number of supply and demand offers},
symbolic x coords={100, 250, 500, 1000},
xtick=data,
bar width=10pt,
scatter/position=absolute,
        node near coords style={
            at={(axis cs:\pgfkeysvalueof{/data point/x},\pgfkeysvalueof{/pgfplots/ymin})},
            anchor=north,
            yshift={-\pgfkeysvalueof{/pgfplots/major tick length}},
            anchor=west,
            rotate=-90,
            color = black,
        },
]
\addplot[blue, fill] coordinates {(100, 54.8) (250, 95.2) (500, 99.2) (1000, 100)};
\end{axis}
\end{tikzpicture}
    \caption{Optimality of the aggregation based method with $\varepsilon = 0.1\%$ threshold over 5 time periods. The number of block bids was 10\% of the number of supply-demand offers.}
    \label{results1}
\end{figure}

\begin{figure}[h!]
    \centering
\begin{tikzpicture}
\pgfplotsset{height=7cm, width=\textwidth,compat=1.17}
\begin{axis}[
ybar,
enlarge y limits={0.2},
enlarge x limits={0.5},
legend style={at={(0.5,1), anchor=north west},
anchor=north,legend columns=-1},
ylabel={Hit rate [\%]},
xlabel={Number of supply and demand offers},
symbolic x coords={100, 250, 500, 1000},
xtick=data,
bar width=10pt,
scatter/position=absolute,
        node near coords style={
            at={(axis cs:\pgfkeysvalueof{/data point/x},\pgfkeysvalueof{/pgfplots/ymin})},
            anchor=north,
            yshift={-\pgfkeysvalueof{/pgfplots/major tick length}},
            anchor=west,
            rotate=-90,
            color = black,
        },
]
\addplot[blue, fill] coordinates {(100, 53.6) (250, 88.4) (500, 100) (1000, 100)};
\end{axis}
\end{tikzpicture}
    \caption{Optimality of the aggregation based method with $\varepsilon = 0.1\%$ threshold over 5 time periods. The number of block bids was 20\% of the number of supply-demand offers.}
    \label{results2}
\end{figure}

According to several tests, as the number of bids increases, the success rates also show an increasing trend (reaches nearly optimum) and the average relative deviation of TSW values decrease in the case of the nominal aggregation based method. This can be explained by the fact that for multiple bids, the clustering algorithm that determines the original bids to be aggregated into a single bid result in larger clusters. When these larger clusters are decomposed, more bids are available for the algorithm in the second step to optimize, thus the chance is lower to reach sub-optimal result.

\subsection{Computational improvement}
Figures \ref{results3} and \ref{results4} depict the computational time for different market configurations.
In real markets, it is usual to have thousands of offers submitted for a single hour (see e.g. \cite{HUPX}).

\begin{figure}[h!]
\centering
\begin{tikzpicture}
\pgfplotsset{height=7cm, width=\textwidth,compat=1.17}
\begin{axis}[
ybar,
enlarge y limits={0.2},
enlarge x limits={0.5},
legend style={at={(0.5,1), anchor=north west},
anchor=north,legend columns=-1},
ylabel={Computational time [s]},
xlabel={Number of supply and demand offers},
symbolic x coords={1000, 5000, 10000},
xtick=data,
bar width=10pt,
scatter/position=absolute,
        node near coords style={
            at={(axis cs:\pgfkeysvalueof{/data point/x},\pgfkeysvalueof{/pgfplots/ymin})},
            anchor=north,
            yshift={-\pgfkeysvalueof{/pgfplots/major tick length}},
            anchor=west,
            rotate=-90,
            color = black,
        },
]
\addplot[blue, fill] coordinates {(1000, 13.4) (5000,162.6) (10000,673)};
\addplot[red, fill] coordinates {(1000, 1.99) (5000, 14.91) (10000, 58.76)};

\legend{Without aggregation, With aggregation}
\end{axis}
\end{tikzpicture}
\caption{Computational time for different bid numbers and solution methods over 10 time periods. The number of block bids was 5\% of the number of supply-demand offers.}
\label{results3}
\end{figure}
 The aggregation based method is faster than the original method by 85\%, 91\%, 91\% for the increasing bid numbers, respectively.

\begin{figure}[h!]
\centering
\begin{tikzpicture}
\pgfplotsset{height=7cm, width=\textwidth,compat=1.17}
\begin{axis}[
ybar,
enlarge y limits={0.2},
enlarge x limits={0.5},
legend style={at={(0.5,1), anchor=north west},
anchor=north,legend columns=-1},
ylabel={Computational time [s]},
xlabel={Number of supply and demand offers},
symbolic x coords={1000, 5000, 10000},
xtick=data,
bar width=10pt,
scatter/position=absolute,
        node near coords style={
            at={(axis cs:\pgfkeysvalueof{/data point/x},\pgfkeysvalueof{/pgfplots/ymin})},
            anchor=north,
            yshift={-\pgfkeysvalueof{/pgfplots/major tick length}},
            anchor=west,
            rotate=-90,
            color = black,
        },
]
\addplot[blue, fill] coordinates {(1000, 17.76) (5000,150.8) (10000,775.6)};
\addplot[red, fill] coordinates {(1000, 1.92) (5000, 15.77) (10000, 57.96)};

\legend{Without aggregation, With aggregation}
\end{axis}
\end{tikzpicture}
\caption{Computational time for different bid numbers and solution methods over 10 time periods. The number of block bids was 10\% of the number of supply-demand offers.}
\label{results4}
\end{figure}
 The aggregation based method is faster than the original method by 89\%, 90\%, 92\% for the increasing bid numbers, respectively.

It is noticeable that with the increase in the number of bids, the aggregation based method outperforms the the standard method used as reference in terms of time requirement.

The figures also show that the amount of block bids slightly influences the optimization time of the nominally aggregated bids. This can be explained by the fact that block bids are not subject to aggregation. In addition, block bids induce binary variables in the problem, which increase the computational time.

\section{Conclusions and future work}
In this work we proposed an approximation-based approach for the clearing of electricity markets. The principle of the approach is the concept of bid-aggregation.
We analyzed the computational gain and optimality of the method in different scenarios regarding the number of submitted bids and the ratio of block bids among them.
Regarding the computational gain, we found that in the case of large bid quantities the method reduced the computational time by one order of magnitude.
Corresponding to the optimality/suboptimality of the proposed method, the results showed that if we assume a large number of bids, the error of the approximation is practically always below 0.5 percent, and if the number of bids was at least 1000-1000 (supply and demand), the error of approximation was below 0.1\% in 99\% of the simulated cases.

\subsection{Future work}
In the future we plan to carry out more simulations to evaluate further details of the method. In addition, the current model is TSW-oriented, and does not explicitly consider market clearing prices (MCPs). Introducing MCPs in the model formulation implies implication type constraints connecting MCPs and bid acceptance constraints, which may be implemented by the big-M method, implying auxiliary binary variables (see e.g. \cite{csercsik2019uncertain}). We expect that the proposed method may be applied efficiently in this case as well, based on the consideration that the significant reduction of bid number leads to the reduction of the number of these auxiliary binary variables.

As mentioned earlier, clearing based on the aggregation of bids can potentially reach a sub-optimal solution. To overcome this issue the algorithm may be initiated with several possible aggregation configurations and the version with the highest TSW value may be considered as the final solution. The optimization problems resulting from these different initial aggregation scenarios may be evaluated in a parallel manner.

\section{Acknowledgements}
{This work has been supported by the \textit{Integrated researcher training program in the disciplines of information technology and computer science} (EFOP-3.6.3-VEKOP-16-2017-00002) project, and
 by the Funds PD 123900 and K 131545 of the Hungarian National Research, Development and Innovation Office. D\{a}vid Csercsik is the grantee of the J\'{a}nos Bolyai Research Scholarship of the Hungarian Academy of Sciences.

\section*{Bibliography}
 \bibliographystyle{elsarticle-num}