# Convergence rates for ensemble-based solutions to optimal control of
uncertain dynamical systems

Olena Melnikov11footnotemark: 1 Johannes Milz H. Milton Stewart School of
Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta,
Georgia 30332<EMAIL_ADDRESS><EMAIL_ADDRESS>
(July 25, 2024)

###### Abstract

We consider optimal control problems involving nonlinear ordinary differential
equations with uncertain inputs. Using the sample average approximation, we
obtain optimal control problems with ensembles of deterministic dynamical
systems. Leveraging techniques for metric entropy bounds, we derive non-
asymptotic Monte Carlo-type convergence rates for the ensemble-based
solutions. Our theoretical framework is validated through numerical
simulations on a harmonic oscillator problem and a vaccination scheduling
problem for epidemic control under model parameter uncertainty.

###### keywords:

stochastic optimization, sample average approximation, Monte Carlo sampling,
sample complexity, uncertain dynamical systems, optimal control

AMS subject classifications. 90C15, 90C30, 49K45, 49N60

##  1 Introduction

We consider the optimal control problem

$\displaystyle\min_{u\in
L^{2}(0,1;\mathbb{R}^{m})}\,\mathbb{E}[F(x^{u}(1,\xi),\xi)]+\psi(u),$ (1)

where for each parameter $\xi\in\Xi$ and control $u(\cdot)\in
L^{2}(0,1;\mathbb{R}^{m})$, $x^{u}(\cdot,\xi)=x(\cdot,\xi)$ solves the
parameterized affine-control dynamical system

$\displaystyle\dot{x}(t,\xi)=f_{0}(x(t,\xi),\xi)+f_{1}(x(t,\xi),\xi)u(t)\quad\text{for
a.e. }t\in(0,1),\quad x(0,\xi)=x_{0}(\xi),$ (2)

where $\Xi$ is a complete separable metric space equipped with its Borel
sigma-algebra, $F\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}$,
$f_{0}\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}^{n}$ and
$f_{1}\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}^{n\times m}$ are Carathéodory
mappings, and $x_{0}\colon\Xi\to\mathbb{R}^{n}$ is measurable. The function
$\psi\colon L^{2}(0,1;\mathbb{R}^{m})\to(-\infty,\infty]$ is proper and
strongly convex with parameter $\alpha>0$. The parameterized initial value
problem in (2) models for uncertain right-hand sides and initial values. With
some abuse notation, we use $\xi$ to denote elements of $\Xi$ and a random
element taking values in $\Xi$.

We use the sample average approximation (SAA) approach [8] to approximate the
infinite dimensional optimization problem (1). Throughout the text, let
$\xi^{1},\xi^{2},\ldots$ be independent and identically distributed
$\Xi$-valued random elements defined on a complete probability space such that
each $\xi^{i}$ has the same distribution as $\xi$. We obtain the SAA problem

$\displaystyle\min_{u\in
L^{2}(0,1;\mathbb{R}^{m})}\,\frac{1}{N}\sum_{i=1}^{N}F(x^{u}(1,\xi^{i}),\xi^{i})+\psi(u).$
(3)

The optimization problem (3) is an optimal control problem with an ensemble of
$N$ dynamical systems. For the numerical solution of (3), we generate a
realization of $\xi^{1},\xi^{2},\ldots,\xi^{N}$ and solve (3) for these fixed
elements.

###  1.1 Contributions

We summarize our main contributions. Our summary uses standard notation, which
is formally introduced in section 2. We define parameterized integrand

$\displaystyle T(u,\xi)\coloneqq F(x^{u}(1,\xi),\xi).$

Furthermore, we define
$\psi_{\alpha}(u)\coloneqq\psi(u)-(\alpha/2)\|u\|_{L^{2}(0,1;\mathbb{R}^{m})}^{2}$,

$\displaystyle
g(u)\coloneqq\mathbb{E}[T(u,\xi)],\quad\text{and}\quad\widehat{g}_{N}(u)\coloneqq\frac{1}{N}\sum_{i=1}^{N}T(u,\xi^{i}).$
(4)

  1. (i)

We establish nonasymptotic mean convergence rates for the SAA optimal values.
Specifically, we show that for all $N\in\mathbb{N}$,

$\displaystyle\mathbb{E}[|\widehat{v}_{N}^{*}-v^{*}|]\leq\frac{\text{Const}}{\sqrt{N}}\bigg{(}1+\frac{1}{\sqrt{\alpha}}\bigg{)},$
(5)

where $v^{*}$ is the optimal value of (1) and $\widehat{v}_{N}^{*}$ is that of
(3). Moreover, Const is a constant that does not depend on the sample size $N$
nor on the strong convexity parameter $\alpha$. However, it can depend on
other problem data, such as the control’s dimension $m$.

  2. (ii)

We demonstrate nonasymptotic mean convergence rates for a criticality measure
for (1) evaluated at SAA critical points: for each critical point
$u_{N}^{*}\in\mathrm{dom}(\psi)$ of (3), that is,
$\widehat{\chi}_{N}(u_{N}^{*})=0$, we show that for all $N\in\mathbb{N}$,

$\displaystyle\mathbb{E}[\chi(u_{N}^{*})]\leq\frac{\text{Const}}{\sqrt{N}}\bigg{(}1+\frac{1}{\sqrt{\alpha}}\bigg{)},$
(6)

where $\mathrm{Const}$ is as in (5) and the criticality measures $\chi$ and
$\widehat{\chi}_{N}$ are defined by

$\displaystyle\chi(u)\coloneqq\|u-\mathrm{prox}_{\psi_{\alpha}}(u-\nabla
g(u)-\alpha u)\|_{L^{2}(0,1;\mathbb{R}^{m})}$ (7)

and

$\displaystyle\widehat{\chi}_{N}(u)\coloneqq\|u-\mathrm{prox}_{\psi_{\alpha}}(u-\nabla\widehat{g}_{N}(u)-\alpha
u)\|_{L^{2}(0,1;\mathbb{R}^{m})}.$

  3. (iii)

Our contributions in (i) and (ii) are derived using uniform expectation bounds
for Hilbert space-valued “sub-Gaussian-type” Carathéodory mappings. Since this
result may be of independent interest, we summarize it here. Let $H$ be a
real, separable Hilbert space, let $C\subset H$ be a nonempty, closed set with
diameter $D>0$, and let $G\colon C\times\Xi\to H$ be a Carathéodory mapping
with $\mathbb{E}[G(x,\xi)]=0$ for each $x\in C$. Suppose that there exists a
constant $M>0$ such that the sub-Gaussian-type estimate

$\displaystyle\mathbb{E}[\cosh(\lambda\|G(x,\xi)-G(y,\xi)\|_{H})]\leq\exp(M^{2}\lambda^{2}\|x-y\|_{H}^{2}/2)\;\text{for
all}\;x,y\in C,\;\lambda\geq 0$

holds. Under these technical conditions, we establish the uniform expectation
bound

$\displaystyle\mathbb{E}[\sup_{x\in
C}\,\|\widehat{G}_{N}(x)\|_{H}]\leq\mathbb{E}[\|\widehat{G}_{N}(x_{0})\|_{H}]+\frac{4\sqrt{3}M}{\sqrt{N}}\int_{0}^{D/2}\sqrt{\ln(2\mathcal{N}(C,\varepsilon))}\,\mathrm{d}\varepsilon,$
(8)

where $\widehat{G}_{N}(x)\coloneqq(1/N)\sum_{i=1}^{N}G(x,\xi^{i})$, and
$\mathcal{N}(C,\varepsilon)$ is the $\varepsilon$-covering number of $C$ (see
section 2 for a definition). This uniform expectation bound is valid for each
fixed $x_{0}\in C$.

  4. (iv)

Contributions (i) through (iii) are technical in nature. Our final
contribution involves formulating vaccination scheduling for epidemic control
with random inputs as a risk-neutral optimal control problem. We validate our
theoretical findings through numerical simulations for this vaccination
scheduling problem.

###  1.2 Related work

The SAA approach has been used for various optimal control problems with
uncertain dynamics, notably in optimizing batch reactors with model parameter
uncertainty [24, 26]. Recent texts have established some of its theoretical
properties, including asymptotic consistency and sample complexity. The paper
[19] analyzes the almost sure epiconvergence of the SAA objective function and
the weak consistency of SAA critical points. Further asymptotic properties are
demonstrated in [25]. Moreover, the work [14] establishes sample complexity of
strongly convex stochastic optimization in Hilbert space. These contributions
broaden the theoretical justification of the SAA approach, extending its
utility to a wider range of infinite dimensional decision-making tasks notably
those constrained by dynamical systems.

Additionally, theoretical analyses of the SAA framework has been extended to
other infinite dimensional stochastic nonconvex optimization problems with
differential equation constraints. This includes the derivation of sample size
estimates for semilinear PDE-constrained optimization under uncertainty [18]
and the consistency analysis of sample-based approximations to risk-averse
stochastic optimization problem with infinite dimensional decision spaces
[17]. For risk-neutral linear elliptic PDE-constrained optimization, mean
convergence rates for SAA optimal values and solutions, as well as a central
limit theorem, are established in [22]. Parts of our theoretical analysis is
based on the construction of a deterministic compact set containing the SAA
critical points. This construction is inspired by the analysis developed in
[18], but we use different statistical tools to derive convergence rates.

The right-hand side of (8) resembles the classical metric entropy bound for
the suprema of sub-Gaussian processes (see, e.g., [5, Thm. 3.1 on p. 95]).
Metric entropy bounds are one of the tools for analyzing the SAA approach;
see, e.g., [12].

###  1.3 Manuscript overview

We present notation and collect basic terminology in section 2. In section 3,
we formulate technical assumptions on the optimal control problem and the
initial value problem. These assumptions enable us to establish the gradient
and control regularity in sections 4 and 5. Our main results on the
convergence rates for the SAA optimal values and critical points are presented
in sections 6 and 7. These theoretical results are empirically validated in
section 8 through two examples. Our appendices collect and establish results
necessary for establishing the SAA convergence rates. While appendix A states
essentially known first-order optimality conditions, appendix B derives bounds
on the covering numbers of square-integrable mappings with square-integrable
derivatives. Finally, appendix C concludes by establishing the uniform
exponential tail bounds.

##  2 Notation

We collect basic notation and terminology. If not stated otherwise, equations
involving random elements hold with probability one. Let $H$ be a Hilbert
space with inner product $(\cdot,\cdot)_{H}$. We use $\|\cdot\|_{H}$ to denote
the norm induced by its inner product $(\cdot,\cdot)_{H}$. For a Banach space
$X$, $\|\cdot\|_{X}$ denotes its norm. For Banach spaces $X_{1}$, $\ldots$,
$X_{m}$, we equip the Cartesian product $X_{1}\times\cdots\times X_{m}$ with
the norm $\|(x_{1},\ldots,x_{m})\|_{X_{1}\times\cdots\times
X_{m}}\coloneqq(\sum_{i=1}^{m}\|x_{i}\|_{X_{i}}^{2})^{1/2}$ if $X_{1}$,
$\ldots$, $X_{m}$ are Hilbert spaces, and
$\|(x_{1},\ldots,x_{m})\|_{X_{1}\times\cdots\times
X_{m}}\coloneqq\sum_{i=1}^{m}\|x_{i}\|_{X_{i}}$ otherwise. We denote by
$\mathbb{B}_{X}(0;r)$ the open unit norm ball with radius $r$ and by
$\overline{\mathbb{B}}_{X}(0;r)$ its closure. If not stated otherwise, we
equip $\mathbb{R}^{m}$ with the Euclidean norm $\|\cdot\|_{2}$. We denote the
spectral norm of a matrix $A$ by $\|A\|_{2}$. We call a function $G$ defined
on a Hilbert space $H$ strongly convex with parameter $\alpha>0$ if
$G(\cdot)-(\alpha/2)\|\cdot\|_{H}^{2}$ is convex. The domain of a function
$\phi$ is denoted by $\mathrm{dom}(\phi)$. For a proper convex function $\phi$
on a Hilbert space $H$, its proximity operator is defined by
$\mathrm{prox}_{\phi}(x)\coloneqq\mathrm{argmin}_{y\in
H}\,(1/2)\|x-y\|_{H}^{2}+\phi(y)$.

The spaces $L^{p}(0,1;\mathbb{R}^{m})$ $(1\leq p\leq\infty)$ denote the
“usual” Lebesgue spaces equipped with their canonical norms. For $1\leq
p\leq\infty$, we define the Sobolev–Bochner space

$\displaystyle W_{1}^{p}(0,1;\mathbb{R}^{m})\coloneqq\\{\,v\in
L^{p}(0,1;\mathbb{R}^{m})\colon\,\dot{v}\in L^{p}(0,1;\mathbb{R}^{m})\,\\}.$

Here, $\dot{v}$ denotes the distributional (time) derivative of $v$. We equip
$W_{1}^{p}(0,1;\mathbb{R}^{m})$ with
$\|v\|_{W_{1}^{p}(0,1;\mathbb{R}^{m})}\coloneqq\big{(}\|v\|_{L^{p}(0,1;\mathbb{R}^{m})}^{p}+\|\dot{v}\|_{L^{p}(0,1;\mathbb{R}^{m})}^{p}\big{)}^{1/p}$
if $1\leq p<\infty$ and
$\|v\|_{W_{1}^{\infty}(0,1;\mathbb{R}^{m})}\coloneqq\|v\|_{L^{\infty}(0,1;\mathbb{R}^{m})}+\|\dot{v}\|_{L^{\infty}(0,1;\mathbb{R}^{m})}$
if $p=\infty$.

Let $H_{2}\subset H_{1}$ be convex sets in a real Hilbert space $H$. Adapting
Definition 5.1.26 in [21], we say that $G\colon H_{1}\to\mathbb{R}$ is Gateaux
differentiable relative to $H_{1}$ at a point $x\in H_{2}$ if there exists a
linear bounded operator $G^{\prime}(x)$ on $H$ such that for each $h\in H$
with $x+\lambda h\in H_{1}$ for some $\lambda>0$, $G$ is directionally
differentiable at $x$ in the direction $h$ with directional derivative
$G^{\prime}(x)h$. We say that $G$ is Gateaux differentiable on $H_{2}$
relative to $H_{1}$ if it is Gateaux differentiable relative to $H_{1}$ at
each $x\in H_{2}$. We denote by $\nabla G(x)$ the Riesz representation of
$G^{\prime}(x)$. (Different from Definition 5.1.26 in [21], we do not require
$G$ be continuous). We denote the Gateaux derivative of a function $G$ by
$DG$, $D_{x}G$, or $G^{\prime}$.

For a metric space $(Y,d_{Y})$ and a nonempty, totally bounded subset
$Y_{0}\subset Y$, we define the $\nu$-covering number
$\mathcal{N}(Y_{0},Y,\nu)$ as the minimal number of points in a $\nu$-net of
$Y_{0}$ in $Y$. We note that the elements in $\nu$-nets of $Y_{0}$ in $Y$ are
not required to be contained in $Y_{0}$. Since $Y_{0}$ is totally bounded,
$(Y_{0},d_{Y})$ is a totally bounded metric space. For abbreviation, we define
$\mathcal{N}(Y_{0},\nu)\coloneqq\mathcal{N}(Y_{0},Y_{0};\nu)$.

##  3 Assumptions on optimal control problem

We state a number of assumptions on the optimal control problem (1) and the
initial value problem (2). The assumptions ensure, e.g, well-posedness of the
optimal control problems, differentiability of sample average and expectation
functions, and existence and uniqueness of solutions to the parameterized
initial value problem. Moreover, these conditions allow us to demonstrate that
the integrand in (1) has square integrable time derivatives for each parameter
$\xi$. These conditions are mainly inspired by those used in [19, 21].

We define

$\displaystyle
r_{\psi}\coloneqq\sup_{u\in\mathrm{dom}(\psi)}\,\|u\|_{L^{\infty}(0,1;\mathbb{R}^{m})}.$

Our analysis is based on the following conditions on the control regularizer.

######  Assumption 1.

The function $\psi\colon L^{2}(0,1;\mathbb{R}^{m})\to(-\infty,\infty]$ is
strongly convex with parameter $\alpha>0$. Moreover, $1\leq r_{\psi}<\infty$.

Assumption 1 ensures that the domain of $\psi$ is a bounded subset in
$L^{\infty}(0,1;\mathbb{R}^{m})$. Assumption 1 is, e.g., satisfied for the
sparsity-promoting control regularizer

$\psi(u)=(\alpha/2)\|u\|_{L^{2}(0,1;\mathbb{R}^{m})}^{2}+\beta\|u\|_{L^{1}(0,1;\mathbb{R}^{m})}+I_{\\{u\in
L^{2}(0,1;\mathbb{R}^{m})\colon\|u\|_{L^{\infty}(0,1;\mathbb{R}^{m})}\leq
1\\}}(u),$

where $\alpha>0$, and $\beta\geq 0$. Moreover, $I_{Y}\colon
L^{2}(0,1;\mathbb{R}^{m})\to\\{0,\infty\\}$ is the indicator function of
$Y\subset L^{2}(0,1;\mathbb{R}^{m})$, that is, $I(u)=0$ if $u\in Y$ and
$I(u)=\infty$ otherwise.

The next assumption is adapted from Assumption 1 in [19]. It requires that the
state be contained in a deterministic compact set.

######  Assumption 2.

The set $X_{0}\subset\mathbb{R}^{n}$ with $0\in X_{0}$ is convex and compact
such that for each $\xi\in\Xi$ and
$u\in\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})$, it holds that
$x^{u}(t,\xi)\in X_{0}$ for a.e. $t\in[0,1]$.

Assumption 2 implies that the parameterized initial value is (essentially)
bounded. More precisely, Assumption 2 ensures the existence of a constant
$V_{0}\in[1,\infty)$ such that $\|x_{0}(\xi)\|_{2}\leq V_{0}$ for all
$\xi\in\Xi$.

The next assumption is motivated by Assumption 5.6.2 in [21] and Assumptions
2, 3, and 5 in [19]. We define the dynamical system’s right-hand side
$f\colon\mathbb{R}^{n}\times\mathbb{R}^{m}\times\Xi\to\mathbb{R}^{n}$ by

$\displaystyle f(x,u,\xi)\coloneqq f_{0}(x,\xi)+f_{1}(x,\xi)u.$ (9)

######  Assumption 3.

Let $r_{\psi}$ be as in Assumption 1, and let $X_{0}$ be as in Assumption 2.

  1. (a)

The mappings $f_{0}\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}^{n}$ and
$f_{1}\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}^{n\times m}$ are Carathéodory
mappings.

  2. (b)

For each $\xi\in\Xi$, $f_{0}(\cdot,\xi)$ and $f_{1}(\cdot,\xi)$ are
continuously differentiable.

  3. (c)

There exists a constant $L_{f}^{\prime}\in[1,\infty)$ such that
$f(\cdot,\cdot,\xi)$, $f_{1}(\cdot,\xi)$, and $D_{x}f(\cdot,\cdot,\xi)$ are
Lipschitz continuous on an open neighborhood of
$X_{0}\times\\{u\in\mathbb{R}^{m}\colon\|u\|_{2}\leq 2r_{\psi}\\}$ with
Lipschitz constant $L_{f}^{\prime}$ for each $\xi\in\Xi$. Moreover,
$\|f(0,0,\xi)\|_{2}\leq L_{f}^{\prime}$ for all $\xi\in\Xi$.

  4. (d)

The function $F\colon\mathbb{R}^{n}\times\Xi\to\mathbb{R}$ is a Carathéodory
function, and $F(\cdot,\xi)$ is continuously differentiable for each
$\xi\in\Xi$. There exists a constant $L_{F}^{\prime}\in[1,\infty)$ such that
$F(\cdot,\xi)$ and $D_{x}F(\cdot,\xi)$ are Lipschitz continuous on an open
neighborhood of $X_{0}$ with Lipschitz constant $L_{F}^{\prime}$ for each
$\xi\in\Xi$.

  5. (e)

The parameterized initial value $x_{0}\colon\Xi\to\mathbb{R}^{n}$ is
measurable.

Assumptions 1, 3 and 2 ensure for all $x\in X_{0}$, $\xi\in\Xi$, and
$u\in\mathbb{R}^{m}$ with $\|u\|_{2}\leq r_{\psi}$,

$\displaystyle\begin{aligned} \|f(x,u,\xi)\|_{2}&\leq
2L_{f}^{\prime}r_{\psi}(\|x\|_{2}+1),\quad\|D_{x}f(x,u,\xi)\|_{2}\leq
L_{f}^{\prime},\quad\|f_{1}(x,\xi)\|_{2}\leq L_{f}^{\prime},\\\
\|D_{x}f_{1}(x,\xi)\|_{2}&\leq L_{f}^{\prime},\end{aligned}$ (10)

where the third inequality is a consequence of $D_{u}f(x,u,\xi)=f_{1}(x,\xi)$
and the Lipschitz continuity of $f(\cdot,\cdot,\xi)$.

##  4 Gradient regularity via adjoint and state regularity

The section’s main purpose is to show that the objective function’s gradient
is contained in a compact set of the feasible set. The analysis performed is
based on the results established in [19, 21].

We recall that the control problem’s integrand
$T\colon\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})\times\Xi\to[0,\infty)$
is defined by

$\displaystyle T(u,\xi)\coloneqq F(x^{u}(1,\xi),\xi).$

If Assumptions 1, 3 and 2 hold and $\xi\in\Xi$, then Corollary 5.6.9 in [21]
ensures that $T(\cdot,\xi)$ is Gateaux differentiable on $\mathrm{dom}(\psi)$
relative to $\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})$ with

$\displaystyle[\nabla_{u}T(u,\xi)](t)=f_{1}(x^{u}(t,\xi),\xi)^{T}p^{u}(t,\xi)\quad\text{for
a.e. }t\in(0,1),$ (11)

where $p^{u}(\cdot,\xi)=p(\cdot,\xi)$ solves for a.e. $t\in(0,1)$ the
parameterized adjoint equation

$\displaystyle\dot{p}(t,\xi)=-D_{x}f(x^{u}(t,\xi),u(t),\xi)^{T}p(t,\xi),\quad
p(1,\xi)=\nabla_{x}F(x^{u}(1,\xi),\xi).$ (12)

As mentioned above, the section’s main purpose is to show that
$\nabla_{u}T(u,\xi)$ is contained in a ball in $W_{1}^{2}(0,1;\mathbb{R}^{m})$
for each $(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$. This result is formulated in
the following theorem, with an explicit expression for the ball’s radius.

######  Theorem 4.1.

If Assumptions 1, 3 and 2 hold, then for each
$(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$, we have

$\displaystyle\nabla_{u}T(u,\xi)\in\overline{\mathbb{B}}_{W_{1}^{\infty}(0,1;\mathbb{R}^{m})}(0;R)\subset\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;R),$

where

$\displaystyle R\coloneqq
R_{0}+6(L_{f}^{\prime})^{2}(1+V_{0})\mathrm{e}^{2L_{f}^{\prime}r_{\psi}}R_{0}\quad\text{with}\quad
R_{0}$ $\displaystyle\coloneqq
L_{f}^{\prime}(1+L_{F}^{\prime})\mathrm{e}^{L_{f}^{\prime}}.$ (13)

We prepare our proof of Theorem 4.1. The following lemma collects basic
properties $T$ mainly from [19]. We define

$\displaystyle
L_{T}^{\prime}\coloneqq\sqrt{2}L_{F}^{\prime}L_{f}^{\prime}\mathrm{e}^{L_{f}^{\prime}}.$
(14)

######  Lemma 4.2.

Let Assumptions 1, 3 and 2 hold. Then the following statements are true.

  1. (a)

For each $\xi\in\Xi$, $T(\cdot,\xi)$ is Lipschitz continuous with Lipschitz
constant $L_{T}^{\prime}$, and Gateaux differentiable on $\mathrm{dom}(\psi)$
relative to $\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})$.

  2. (b)

There exists a constant $L_{\nabla T}^{\prime}\in[1,\infty)$ such that
$\nabla_{u}T(\cdot,\xi)$ is Lipschitz continuous on $\mathrm{dom}(\psi)$ with
Lipschitz constant $L_{\nabla T}^{\prime}$ for each $\xi\in\Xi$.

  3. (c)

The mappings $T$ and $\nabla_{u}T$ are Carathéodory mappings on
$\mathrm{dom}(\psi)\times\Xi$.

  4. (d)

The function $u\mapsto\mathbb{E}[T(u,\xi)]$ is Gateaux differentiable on
$\mathrm{dom}(\psi)$ relative to
$\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})$.

###### Proof.

  1. (a)

The Lipschitz statement follows from an application of Lemma 5.6.3 in [21].
See also the proof of Lemma 5.6.7 in [21], and Lemma 3.5 in [19]. Furthermore,
Corollary 5.6.9 in [21] ensures the Gateaux differentiability statement.

  2. (b)

This statement follows from an application of Lemma 5.2 in [19].

  3. (c)

According to Lemma 3.5 in [19], $T$ is a Carathéodory mapping. Using Lemma 2.3
in [19], composition rules, and (11), we can show that $\nabla_{u}T(u,\cdot)$
on $\Xi$ is measurable. See also Lemma 5.2 in [19].

  4. (d)

An application of Lemma 4.3 in [19] yields the assertion.

∎

Using chain and product rules (cf. [11, p. 340] and [1, Cor. 3.2]), and (11),
we find that for a.e. $t\in(0,1)$,

$\displaystyle[\nabla_{u}T(u,\xi)]^{\prime}(t)=D_{t}[f_{1}(x^{u}(t,\xi),\xi)^{T}]p^{u}(t,\xi)+f_{1}(x^{u}(t,\xi),\xi)^{T}\dot{p}^{u}(t,\xi).$
(15)

This identity allows us to establish stability estimates on the gradient’s
time derivative.

######  Lemma 4.3.

Let Assumptions 1, 3 and 2 hold. Then the following statements are true.

  1. (a)

For each
$(u,\xi)\in\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})\times\Xi$,
the dynamical system (2) has a unique solution $x^{u}(\cdot,\xi)\in
W_{1}^{\infty}(0,1;\mathbb{R}^{n})$, and for a.e. $t\in[0,1]$,

$\displaystyle\begin{aligned}
\|x^{u}(t,\xi)\|_{2}&\leq(1+V_{0})\mathrm{e}^{2L_{f}^{\prime}r_{\psi}},\\\
\|\dot{x}^{u}(t,\xi)\|_{2}&\leq
2L_{f}^{\prime}r_{\psi}\big{(}\|x^{u}(t,\xi)\|_{2}+1\big{)}.\end{aligned}$
(16)

  2. (b)

For each
$(u,\xi)\in\mathbb{B}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;2r_{\psi})\times\Xi$,
the adjoint system (12) has a unique solution $p^{u}(\cdot,\xi)\in
W_{1}^{\infty}(0,1;\mathbb{R}^{n})$, and for a.e. $t\in[0,1]$,

$\displaystyle\begin{aligned}
\|p^{u}(t,\xi)\|_{2}&\leq(1+L_{F}^{\prime})\mathrm{e}^{L_{f}^{\prime}},\\\
\|\dot{p}^{u}(t,\xi)\|_{2}&\leq
L_{f}^{\prime}\|p^{u}(t,\xi)\|_{2}.\end{aligned}$ (17)

  3. (c)

For all $(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$,

$\displaystyle\|\nabla_{u}T(u,\xi)\|_{L^{\infty}(0,1;\mathbb{R}^{m})}\leq
L_{f}^{\prime}(1+L_{F}^{\prime})\mathrm{e}^{L_{f}^{\prime}}.$

  4. (d)

For all $(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$,

$\displaystyle\|[\nabla_{u}T(u,\xi)]^{\prime}\|_{L^{\infty}(0,1;\mathbb{R}^{m})}\leq
6(L_{f}^{\prime})^{2}(1+L_{F}^{\prime})\mathrm{e}^{L_{f}^{\prime}}(1+V_{0})\mathrm{e}^{2L_{f}^{\prime}r_{\psi}}.$

###### Proof.

  1. (a)

The existence, uniqueness, and the state stability estimate follows from
Proposition 5.6.5 in [21]. Using (2), and (10), we obtain the second stability
estimate.

  2. (b)

The existence and uniqueness of the solution to the adjoint equation follows
from Corollary 5.6.9 in [21].

The mean value theorem ensures $\|D_{x}F(x,\xi)\|_{2}\leq L_{F}^{\prime}$ for
all $(x,\xi)\in X_{0}$. Combined with (10), and Exercise 5.6.6 in [21], we
obtain the first adjoint state stability estimate. The second is a consequence
of the adjoint equation (12).

  3. (c)

Using (10), and (11), we have for all $(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$
and $t\in(0,1)$,

$\displaystyle\|[\nabla_{u}T(u,\xi)](t)\|_{2}\leq\|f_{1}(x^{u}(t,\xi),\xi)\|_{2}\|p^{u}(t,\xi)\|_{2}\leq
L_{f}^{\prime}\|p^{u}(t,\xi)\|_{2}.$

Combined with (17), we obtain the estimate.

  4. (d)

Using (10) and (15), we obtain

$\displaystyle\|[\nabla_{u}T(u,\xi)]^{\prime}\|_{L^{\infty}(0,1;\mathbb{R}^{m})}$
$\displaystyle\leq
L_{f}^{\prime}\|p^{u}(\cdot,\xi)\|_{L^{\infty}(0,1;\mathbb{R}^{m})}\|\dot{x}^{u}(\cdot,\xi)\|_{L^{\infty}(0,1;\mathbb{R}^{m})}$
$\displaystyle\quad+L_{f}^{\prime}\|\dot{p}^{u}(\cdot,\xi)\|_{L^{\infty}(0,1;\mathbb{R}^{m})}.$

Now, (16) and (17) ensure the bound.

∎

The proof of Theorem 4.1 is now a direct consequence of Lemma 4.3.

###### Proof of Theorem 4.1.

Since

$\displaystyle\|\nabla_{u}T(u,\xi)\|_{W_{1}^{\infty}(0,1;\mathbb{R}^{m})}=\|\nabla_{u}T(u,\xi)\|_{L^{\infty}(0,1;\mathbb{R}^{m})}+\|\nabla_{u}T(u,\xi)^{\prime}\|_{L^{\infty}(0,1;\mathbb{R}^{m})}$

Lemma 4.3 ensures the assertion. ∎

##  5 Control regularity via gradient regularity

We show that SAA critical points are contained in a compact deterministic
subset in the control space $L^{2}(0,1;\mathbb{R}^{m})$. The technical results
presented in this section form the basis for establishing our mean convergence
rates in sections 6 and 7. We call $u\in\mathrm{dom}(\psi)$ a critical point
of (1) if $\chi(u)=0$. Moreover, we refer to $u_{N}\in\mathrm{dom}(\psi)$ as a
critical point of the SAA problem if $\widehat{\chi}_{N}(u_{N})=0$.

We define the subset $\mathcal{C}\subset L^{2}(0,1;\mathbb{R}^{m})$ by

$\displaystyle\mathcal{C}\coloneqq\mathrm{prox}_{(1/\alpha)\psi_{\alpha}}(-(1/\alpha)\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;R)).$
(18)

Our next theorem establishes the covering numbers of $\mathcal{C}$. Moreover,
it demonstrates that each critical point of the control problem (1) and of its
SAA problem (3) is contained in $\mathcal{C}$.

######  Theorem 5.1.

If Assumptions 1, 3 and 2 holds, then the following statements are true.

  1. (a)

The diameter of $\mathcal{C}$ in $L^{2}(0,1;\mathbb{R}^{m})$ is at most
$2r_{\psi}$.

  2. (b)

For all $\nu>0$, the $\nu$-covering number of $\mathcal{C}$ with respect to
the $L^{2}(0,1;\mathbb{R}^{m})$-norm satisfies

$\displaystyle\log_{2}(\mathcal{N}(\mathcal{C},\nu))\leq\varrho\sqrt{m}(Rm/(\alpha\nu)),$

where $\varrho>0$ is a constant independent of $m$, $\nu$, $R$, and $\alpha$.

  3. (c)

We have

$\displaystyle\int_{0}^{r_{\psi}}\sqrt{\ln(2\mathcal{N}(\mathcal{C},\varepsilon))}\,\mathrm{d}\varepsilon\leq
2r_{\psi}\sqrt{1+\varrho\sqrt{m}(Rm/\alpha)}.$

  4. (d)

Each critical point of (1) and of (3) is contained in $\mathcal{C}$.

###### Proof.

  1. (a)

We have $\mathcal{C}\subset\mathrm{dom}(\psi_{\alpha})=\mathrm{dom}(\psi)$.
Assumption 1 ensures that
$\mathrm{dom}(\psi)\subset\overline{\mathbb{B}}_{L^{\infty}(0,1;\mathbb{R}^{m})}(0;r_{\psi})$.
Hence the diameter of $\mathcal{C}$ in $L^{2}(0,1;\mathbb{R}^{m})$ is at most
$2r_{\psi}$.

  2. (b)

We define $\varepsilon\coloneqq(\alpha/R)\nu$. Proposition B.1 ensures the
existence of points
$v_{1},\ldots,v_{K}\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;1)$,
such that $\|v-v_{k(v)}\|_{L^{2}(0,1;\mathbb{R}^{m})}\leq\varepsilon$ for all
$v\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;1)$, where
$k(v)\coloneqq\mathrm{argmin}_{k\in\\{1,\ldots,K\\}}\|v-v_{k}\|_{L^{2}(0,1;\mathbb{R}^{m})}$.
Additionally, $K\leq 2^{\varrho\sqrt{m}(m/\varepsilon)}$, with $\varrho>0$
independent of $m$ and $\varepsilon$.

Let $w\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;R)$. Then
$v\coloneqq w/R\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;1)$.
Since proximity operators are nonexpansive,

$\displaystyle\|\mathrm{prox}_{(1/\alpha)\psi_{\alpha}}(-(1/\alpha)w)-\mathrm{prox}_{(1/\alpha)\psi_{\alpha}}(-(1/\alpha)Rv_{k(v)})\|_{L^{2}(0,1;\mathbb{R}^{m})}$
$\displaystyle\quad\leq(1/\alpha)\|w-Rv_{k(v)}\|_{L^{2}(0,1;\mathbb{R}^{m})}\leq(1/\alpha)\|w-Rv_{k(v)}\|_{L^{2}(0,1;\mathbb{R}^{m})}$
$\displaystyle\quad=(R/\alpha)\|v-v_{k(v)}\|_{L^{2}(0,1;\mathbb{R}^{m})}\leq(R/\alpha)\varepsilon=\nu.$

  3. (c)

Combining the previous part, Example C.6, and the bound

$\displaystyle\sqrt{r_{\psi}^{2}+r_{\psi}\varrho\sqrt{m}(Rm/\alpha)}+r_{\psi}\sqrt{\varrho\sqrt{m}(Rm/\alpha)}\leq
2r_{\psi}\sqrt{1+\varrho\sqrt{m}(Rm/\alpha)},$

we obtain the integral bound.

  4. (d)

Let $\bar{u}$ be a critical point of (1). Lemmas A.1 and 4.2 ensure

$\bar{u}=\mathrm{prox}_{(1/\alpha)\psi_{\alpha}}(-(1/\alpha)\nabla\mathbb{E}[T(\bar{u},\xi)]).$

Moreover, Theorem 4.1 and the dominated convergence theorem yield
$\nabla\mathbb{E}[T(\bar{u},\xi)]=\mathbb{E}[\nabla_{u}T(\bar{u},\xi)]\in
L^{2}(0,1;\mathbb{R}^{m})$. Theorem 4.1 also ensures $\nabla_{u}T(u,\xi)\in
W_{1}^{2}(0,1;\mathbb{R}^{m})$ for all
$(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$. Since $W_{1}^{2}(0,1;\mathbb{R}^{m})$
is reflexive, $\overline{\mathbb{B}}_{W_{1}^{2}(0,1;\mathbb{R}^{m})}(0;R)$ is
closed in $L^{2}(0,1;\mathbb{R}^{m})$. Combined with Proposition 1.2.12 in
[7], we find that $\mathbb{E}[\nabla_{u}T(\bar{u},\xi)]\in
W_{1}^{2}(0,1;\mathbb{R}^{m})$. Hence $\bar{u}\in\mathcal{C}$.

Now, let $\bar{u}_{N}$ be a critical point of the SAA problem (3). Lemma A.1
yields
$\bar{u}_{N}=\mathrm{prox}_{(1/\alpha)\psi_{\alpha}}(-(1/\alpha)\nabla\widehat{g}_{N}(\bar{u}_{N}))$.
Since $\nabla_{u}T(u,\xi)\in W_{1}^{2}(0,1;\mathbb{R}^{m})$ for all
$(u,\xi)\in\mathrm{dom}(\psi)\times\Xi$, we have
$\nabla\widehat{g}_{N}(\bar{u}_{N})=(1/N)\sum_{i=1}^{N}\nabla_{u}T(\bar{u}_{N},\xi^{i})\in
W_{1}^{2}(0,1;\mathbb{R}^{m})$. Hence $\bar{u}_{N}\in\mathcal{C}$.

∎

##  6 Convergence rates of optimal values

We establish mean convergence rates for SAA optimal values, our first main
contribution. Throughout the section we assume that the control problems (1)
and (3) have solutions. Specifically, let $u^{*}$ be a solution to (1) and for
each $N\in\mathbb{N}$, let $u_{N}^{*}$ be a solution to (3). Moreover, we
denote by $v^{*}$ the optimal value of (1) and by $\widehat{v}_{N}^{*}$ that
of (3) .

######  Theorem 6.1.

If Assumptions 1, 3 and 2 hold and $u_{0}\in\mathrm{dom}(\psi)$, then for all
$N\in\mathbb{N}$,

$\displaystyle\mathbb{E}[|\widehat{v}_{N}^{*}-v^{*}|]$
$\displaystyle\leq\frac{(\mathbb{E}[(T(u_{0},\xi)-\mathbb{E}[T(u_{0},\xi)])^{2}])^{1/2}}{\sqrt{N}}$
$\displaystyle\quad+\frac{16\sqrt{3}L_{T}^{\prime}r_{\psi}}{\sqrt{N}}\bigg{(}1+\frac{\varrho\sqrt{m}Rm}{\alpha}\bigg{)}^{1/2}.$

###### Proof.

Since each solution to (1) is a critical point of (1), Theorem 5.1 ensures
that the solution set of (1) is contained in $\mathcal{C}$. In particular,
$u^{*}\in\mathcal{C}$. Analogously, each SAA solution is contained in
$\mathcal{C}$. Hence,

$\displaystyle|\widehat{v}_{N}^{*}-v^{*}|\leq\sup_{u\in\mathcal{C}}\bigg{|}\mathbb{E}[T(u,\xi)]-\frac{1}{N}\sum_{i=1}^{N}T(u,\xi^{i})\bigg{|}.$

Next, we estimate the supremum’s expected value using Theorem C.1. For each
$\xi\in\Xi$, Lemma 4.2 ensures that $T(\cdot,\xi)$ is Lipschitz continuous
with Lipschitz constant $L_{T}^{\prime}$ on $\mathrm{dom}(\psi)$. We define
$G(u,\xi)\coloneqq T(u,\xi)-\mathbb{E}[T(u,\xi)]$. For each $\xi\in\Xi$,
$G(\cdot,\xi)$ is Lipschitz continuous with Lipschitz constant
$2L_{T}^{\prime}$ on $\mathrm{dom}(\psi)$. Hence $G$ satisfies (21) with
$M=2L_{T}^{\prime}$. Now Theorems C.1 and C.2 ensure

$\displaystyle\mathbb{E}[|\widehat{v}_{N}^{*}-v^{*}|]$
$\displaystyle\leq\frac{(\mathbb{E}[(T(u_{0},\xi)-\mathbb{E}[T(u_{0},\xi)])^{2}])^{1/2}}{\sqrt{N}}+\frac{8\sqrt{3}L_{T}^{\prime}}{\sqrt{N}}\int_{0}^{r_{\psi}}\sqrt{\ln(2\mathcal{N}(\mathcal{C},\varepsilon))}\,\mathrm{d}\varepsilon.$

Combined with the entropy integral bound provided by Theorem 5.1, we obtain
the assertion. ∎

##  7 Convergence rates for criticality measures

We establish mean convergence rates for SAA critical points. We measure
accuracy of SAA critical points via the criticality measure defined in (7).

######  Theorem 7.1.

Let Assumptions 1, 3 and 2 hold. For each $N\in\mathbb{N}$, let
$u_{N}^{*}\in\mathrm{dom}(\psi)$ be a measurable critical point of the SAA
problem (3). If $u_{0}\in\mathrm{dom}(\psi)$, then for all $N\in\mathbb{N}$,

$\displaystyle\mathbb{E}[\chi(u_{N}^{*})]$
$\displaystyle\leq\frac{(\mathbb{E}\|\nabla_{u}T(u_{0},\xi)-\mathbb{E}[\nabla_{u}T(u_{0},\xi)]\|_{L^{2}(0,1;\mathbb{R}^{m})}^{2})^{1/2}}{\sqrt{N}}$
$\displaystyle\quad+\frac{16\sqrt{3}L_{\nabla
T}^{\prime}r_{\psi}}{\sqrt{N}}\bigg{(}1+\frac{\varrho\sqrt{m}Rm}{\alpha}\bigg{)}^{1/2}.$

###### Proof.

Our proof follows steps similar to those in the proof of Theorem 6.1. For all
$u\in\mathrm{dom}(\psi)$, the nonexpansiveness of proximity operators ensures

$\displaystyle|\chi(u)-\widehat{\chi}_{N}(u)|\leq\|\nabla
g(u)-\nabla\widehat{g}_{N}(u)\|_{L^{2}(0,1;\mathbb{R}^{m})}.$

Combined with $\widehat{\chi}_{N}(u_{N}^{*})=0$ and $u_{N}^{*}\in\mathcal{C}$,
we obtain

$\displaystyle\mathbb{E}[\chi(u_{N}^{*})]\leq\mathbb{E}\big{[}\sup_{u\in\mathcal{C}}\,\|\nabla
g(u)-\nabla\widehat{g}_{N}(u)\|_{L^{2}(0,1;\mathbb{R}^{m})}\big{]}.$

Next, we estimate the above right-hand side using Theorem C.1. For each
$\xi\in\Xi$, Lemma 4.2 ensures that $\nabla_{u}T(\cdot,\xi)$ is Lipschitz
continuous with Lipschitz constant $L_{\nabla T}^{\prime}$ on
$\mathrm{dom}(\psi)$. We define
$G(u,\xi)\coloneqq\nabla_{u}T(u,\xi)-\mathbb{E}[\nabla_{u}T(u,\xi)]$. For each
$\xi\in\Xi$, $G(\cdot,\xi)$ is Lipschitz continuous with Lipschitz constant
$2L_{\nabla T}^{\prime}$ on $\mathrm{dom}(\psi)$. Hence $G$ satisfies (21)
with $M=2L_{\nabla T}^{\prime}$. Combined with Theorems 5.1, C.1, C.1 and C.2,
we obtain the assertion. ∎

##  8 Numerical illustrations

We consider two numerical examples. The main purpose of our numerical
illustrations is to demonstrate the Monte Carlo convergence rate for the SAA
optimal value and SAA critical points established in sections 6 and 7. Before
discussing these examples, we provide details on the discretization of the
control problems, implementation details, and a description of simulation
output. We also introduce some terminology.

We refer to

$\displaystyle\min_{u\in
L^{2}(0,1;\mathbb{R}^{m})}\,F(x^{u}(1,\mathbb{E}[\xi]),\mathbb{E}[\xi])+\psi(u),$

as the nominal problem to (1). It is obtained be replacing the random element
$\xi$ by its expected value and simulating the dynamical system (2) with this
expected value. For both numerical examples, we visualize the nominal solution
and the solution to a certain SAA problem, which is introduced below and
referred to as the reference problem. These graphical illustrations
demonstrate that nominal solutions and reference solutions can have different
characteristics.

To approximate the optimal value of (1) and the criticality measures $\chi$,
we generated $N_{\text{ref}}=4096$ samples using a realization of a scrambled
Sobol’ sequence and transformed these samples to take values in $\Xi$. We
denote them by $\zeta^{1},\ldots,\zeta^{N_{\text{ref}}}$ and formulate the
reference problem

$\displaystyle\min_{u\in
L^{2}(0,1;\mathbb{R}^{m})}\,\frac{1}{N_{\text{ref}}}\sum_{i=1}^{N_{\text{ref}}}F(x^{u}(1,\zeta^{i}),\zeta^{i})+\psi(u).$

Moreover, we approximate the criticality measure $\chi$ by

$\displaystyle\chi_{\text{ref}}(u)\coloneqq\Big{\|}u-\mathrm{prox}_{\psi_{\alpha}}\Big{(}u-\frac{1}{N_{\text{ref}}}\sum_{i=1}^{N_{\text{ref}}}\nabla_{u}T(u,\zeta^{i})-\alpha
u\Big{)}\Big{\|}_{L^{2}(0,1;\mathbb{R}^{m})}.$

The SAA problems’ ensemble of dynamical systems is discretized using multiple
shooting with a fixed step explicit fourth-order Runge–Kutta method. The
control space is discretized using piecewise constant functions defined on a
uniform grid over $[0,t_{f}]$ with $q=50$ subintervals, where $t_{f}>0$ is a
fixed final time. Consequently, the solutions to the SAA problems and the
criticality measures depend on the discretization parameter $q$. For
notational convenience, this discretization parameter is omitted in the
notation. We generated Sobol’ sequences using scipy’s quasi-Monte Carlo
capabilities [23]. To model and solve the discretized SAA problems, we used
the package EnsembleControl [15], which is build on the software framework
CasADi [2], and the optimization solver IPOPT [28]. Our simulation code and
output is archived at [16]. Simulations were performed on a laptop equipped
with a 12th Gen Intel(R) Core(TM) i7-1260P processor and 16 GB of RAM.

###  8.1 Harmonic oscillator

(a)

Figure 1: For the harmonic oscillator problem formulated in section 8.1,
nominal solution _(left)_ and reference solution _(right)_.

(a)  (b)

Figure 2: For the harmonic oscillator problem formulated in section 8.1,
convergence rate of the SAA optimal values _(left)_ and the reference
criticality measure evaluated at SAA critical points _(right)_. The empirical
means $\widehat{\mathbb{E}}$ were computed using $50$ replications. The
convergence rates were computed using least squares.

We consider a harmonic oscillator inspired by the problem formulation in [19,
sect. 7.1]. Our optimization problem is given by

$\displaystyle\min_{u\in
L^{2}(0,1;\mathbb{R}^{2})}\,\mathbb{E}[x_{1}^{u}(1,\xi)^{2}+x_{2}^{u}(1,\xi)^{2}]+\psi(u),$

where for each parameter $\xi\in\Xi\subset\mathbb{R}^{5}$ and control
$u(\cdot)\in L^{2}(0,1;\mathbb{R}^{2})$, $x^{u}(\cdot,\xi)=x(\cdot,\xi)$
solves

$\displaystyle\dot{x}(t,\xi)=\begin{bmatrix}0&-\xi_{1}\\\
\xi_{1}&0\end{bmatrix}x(t,\xi)+\begin{bmatrix}u_{1}(t)\\\
u_{2}(t)\end{bmatrix}+\begin{bmatrix}\xi_{2}\\\ \xi_{3}\end{bmatrix},\quad
t\in(0,1),\quad x(0,\xi)=\begin{bmatrix}1+\xi_{4}\\\ \xi_{5}\end{bmatrix}.$

Moreover, the function $\psi(\cdot)$ is the sum of
$(\alpha/2)\|\cdot\|_{L^{2}(0,1;\mathbb{R}^{2})}^{2}$ and the indicator
function of the set $\\{u\in
L^{2}(0,1;\mathbb{R}^{2})\colon\|u_{i}\|_{L^{\infty}(0,1;\mathbb{R})}\leq
3,\quad i=0,1\\}$. The random variables $\xi_{1},\ldots,\xi_{5}$ are
independent, and uniformly distributed on $[0,2\pi]$, $[-5,5]$, $[-5,5]$,
$[-1/2,1/2]$, and $[-1/2,1/2]$, respectively. The Cartesian product of these
intervals defines the sample space $\Xi$.

Figure 1 depicts the nominal and reference solutions. Figure 2 depicts
convergence rates for the SAA optimal values and critical points. These
empirical rates are close to the theoretical rates established in sections 6
and 7.

###  8.2 Vaccination scheduling for epidemic control under model parameter
uncertainty

(a)

Figure 3: For the vaccination scheduling problem formulated in section 8.2,
nominal solution _(left)_ and reference solution _(right)_.

We formulate an optimal control problem to obtain vaccination schedules for an
epidemic disease under model parameter uncertainty. This problem formulation
is inspired by the optimal control problem formulated in [13, sect. 4] (see
also [10, Chap. 13]). The dynamics is based on an SEIR model described in [13,
sect. 4.2], which tracks the number of individuals who are susceptible,
exposed, infected, and vaccinated, as well as the total population. The model
relies on six key factors: natural death rate, disease death rate, birth rate,
infection incidence, infection rate, and recovery rate. We account for
uncertainty by treating these six parameters as random variables. Our control
variable is the vaccination rate for susceptible individuals. The goal is to
minimize the expected number of infected individuals and the cost of
vaccination.

We consider the optimal control problem

$\displaystyle\min_{u\in
L^{2}(0,t_{f};\mathbb{R})}\,\mathbb{E}\Big{[}\int_{0}^{t_{f}}I^{u}(t,\xi)\mathrm{d}t\Big{]}+\psi(u),$
(19)

where for each control $u(\cdot)\in L^{2}(0,t_{f};\mathbb{R})$ with $0\leq
u(t)\leq 0.9$ for a.e. $t\in(0,1)$ and parameter
$\xi\coloneqq(a,b,c,d,e,g)\in\Xi$, the states $S^{u}(\cdot,\xi)$,
$E^{u}(\cdot,\xi)$, $I^{u}(\cdot,\xi)$, $R^{u}(\cdot,\xi)$, $N^{u}(\cdot,\xi)$
solve the SEIR model (see [13, eqns. (4.1)–(4.5)])

$\displaystyle S^{\prime}(t)$ $\displaystyle=bN(t)-dS(t)-cS(t)I(t)-u(t)S(t),$
$\displaystyle S(0)$ $\displaystyle=S_{0},$ $\displaystyle E^{\prime}(t)$
$\displaystyle=cS(t)I(t)-(e+d)E(t),$ $\displaystyle E(0)$
$\displaystyle=E_{0},$ $\displaystyle I^{\prime}(t)$
$\displaystyle=eE(t)-(g+a+d)I(t),$ $\displaystyle I(0)$ $\displaystyle=I_{0},$
$\displaystyle R^{\prime}(t)$ $\displaystyle=gI(t)-dR(t)+u(t)S(t),$
$\displaystyle R(0)$ $\displaystyle=R_{0},$ $\displaystyle N^{\prime}(t)$
$\displaystyle=(b-d)N(t)-aI(t),$ $\displaystyle N(0)$ $\displaystyle=N_{0}.$

Following [13, p. 76], we choose $\psi(\cdot)$ as the sum of the indicator
function of $\\{u\in L^{2}(0,t_{f};\mathbb{R})\colon 0\leq u(t)\leq 0.9$ a.e.
$t\in(0,1)\\}$ and $(\alpha/2)\|\cdot\|_{L^{2}(0,t_{f};\mathbb{R})}^{2}$ with
$\alpha=2$. For the numerical solution, we removed the fourth state equation,
as none of the other states depend on $R(\cdot)$. The control problem (19) can
be formulated as an instance of (1) via a time transformation and the
introduction of an additional state variable.

We describe the parameter values and our choice of the random vector
$\xi\in\mathbb{R}^{6}$ used for our numerical simulations. As in [13, Table
1], we use the initial states $S_{0}=1000$, $E_{0}=100$, $I_{0}=50$,
$R_{0}=15$. Let $N_{0}=S_{0}+E_{0}+I_{0}+R_{0}$ and let $t_{f}=20$. We choose
the nominal parameter $\bar{\xi}=(0.2,0.525,0.001,0.5,0.5,0.1)$ for
$\xi=(a,b,c,d,e,g)$ using the parameter values in [13, Table 1]. We construct
the random variables $\xi_{j}$, $j=1,\ldots,6$ through random relative
perturbations of the nominal parameter $\bar{\xi}$. Specifically, for
$\sigma=0.15$ and independent random variables $\rho_{j}$, $j=1,\ldots,6$,
each uniformly distributed on the interval $[-1,1]$, we define

$\displaystyle\xi_{j}=(1+\sigma\rho_{j})\bar{\xi}_{j},\quad j=1,\ldots,6.$

Hence, $\Xi$ is chosen as the Cartesian product of the intervals
$[(1-\sigma)\bar{\xi}_{j},(1+\sigma)\bar{\xi}_{j}]$, $j=1,\ldots,6$.

Figure 3 depicts the nominal and reference solutions. Figure 4 depicts
convergence rates for the SAA optimal values and critical points. These
empirical rates closely match the theoretical rates established in sections 6
and 7.

(a)

Figure 4: For the vaccination scheduling problem formulated in section 8.2,
convergence rate of the SAA optimal values _(left)_ and the reference
criticality measure evaluated at SAA critical points _(right)_. The empirical
means $\widehat{\mathbb{E}}$ were computed using $50$ replications. The
convergence rates were computed using least squares.

##  9 Discussion

We have considered risk-neutral optimal control of affine-control dynamics
with random inputs. This manuscript has demonstrated non-asymptotic
convergence rates for SAA optimal values and SAA critical points, which align
with the typical Monte Carlo convergence rate. We have empirically validated
these rates through numerical simulations. An open question remains whether
the dependence on the strong convexity parameter in the convergence rates is
optimal. This consideration may become particularly relevant for applications
of our results to (statistical) inverse problems.

### Acknowledgment

The second author thanks Ashwin Pananjady for insightful discussions on
entropy integral bounds, optimality, and regularization.

### Reproducibility of computational results

Computer code that allows the reader to reproduce the computational results in
this manuscript is available at https://doi.org/10.5281/zenodo.12740932.

##  Appendix A Optimality conditions

We state essentially known first-order necessary optimality conditions for
composite optimization in Hilbert spaces. We use these optimality conditions
to derive first-order optimality conditions for the control problem under
uncertainty and its SAA problem, and to define criticality measures. These
optimality conditions are used to justify the definition of certain
criticality measures.

We consider the composite optimization problem

$\displaystyle\min_{u\in\mathrm{dom}(\psi)}\,G(u)+\psi(u).$ (20)

We formulate assumptions on the optimization problem (20).

######  Assumption 4.

  1. (a)

The space $U$ is a real Hilbert space.

  2. (b)

The function $\psi\colon U\to\mathbb{R}\cup\\{\infty\\}$ is proper, lower
semicontinuous, and strongly convex with parameter $\alpha>0$.

  3. (c)

The set $U_{1}\subset U$ is convex with $\mathrm{dom}(\psi)\subset U_{1}$. The
function $G\colon U_{1}\to\mathbb{R}$ is Gateaux differentiable on
$\mathrm{dom}(\psi)$ relative to $U_{1}$.

Let Assumption 4 hold true. We say that $\bar{u}\in\mathrm{dom}(\psi)$ is a
critical point of (20) if $-\nabla G(\bar{u})\in\partial\psi(\bar{u})$. Here,
$\partial\phi(u)$ denotes the subdifferential of $\phi\colon
U\to(-\infty,\infty]$ at $u\in U$. Let us define
$\psi_{\alpha}(u)\coloneqq\psi(u)-(\alpha/2)\|u\|_{U}^{2}$.

The following essentially known fact characterizes the critical points of
(20).

######  Lemma A.1.

Let Assumption 4 be satisfied. If $\bar{u}\in\mathrm{dom}(\psi)$ is a local
solution to (20), then $-\nabla G(\bar{u})\in\partial\psi(\bar{u})$. If
$\bar{u}\in\mathrm{dom}(\psi)$, then $\bar{u}$ is a critical point of (20) if
and only if

$\bar{u}=\mathrm{prox}_{\gamma\psi_{\alpha}}(\bar{u}-\gamma\nabla
G(\bar{u})-\alpha\gamma\bar{u})\quad\text{for any and all}\quad\gamma>0.$

###### Proof.

The necessary first-order optimality condition is a direct consequence of
well-known statements. Let us define $G_{\alpha}(u)\coloneqq
G(u)+(\alpha/2)\|u\|_{U}^{2}$. We have
$\partial\psi(u)=\partial\psi_{\alpha}(u)+\alpha u$. Therefore, $\bar{u}$ is a
critical point for (20) if and only if $-\nabla
G_{\alpha}(\bar{u})\in\partial\psi_{\alpha}(\bar{u})$ if and only if
$-\gamma\nabla G_{\alpha}(\bar{u})\in\partial[\gamma\psi_{\alpha}(\bar{u})]$.
Combined with [3, Prop. 12.26], we obtain the characterizations of critical
points of (20). ∎

##  Appendix B Covering numbers of vector-valued functions

We provide upper bounds on the covering numbers of
$\overline{\mathbb{B}}_{W_{1}^{2}(0,1,\mathbb{R}^{m})}(0;1)$ with respect to
the $L^{2}(0,1;\mathbb{R}^{m})$-norm. The covering numbers may be nonoptimal
with respect to the space dimension $m$, but are optimal with respect to the
covering radii. In the main text, we use these covering numbers to establish
those of a deterministic set containing all SAA critical points.

######  Proposition B.1.

For some $\varrho>0$ and each $m\in\mathbb{N}$, the binary logarithm of the
$\nu$-covering number of
$\overline{\mathbb{B}}_{W_{1}^{2}(0,1,\mathbb{R}^{m})}(0;1)$ with respect to
the $L^{2}(0,1;\mathbb{R}^{m})$-norm does not exceed $\varrho\sqrt{m}(m/\nu)$
for all sufficiently small $\nu>0$.

###### Proof.

According to Theorem 1.7 in [4], the binary logarithm of the $\nu$-covering
number of $\overline{\mathbb{B}}_{W_{1}^{2}(0,1)}(0;1)$ with respect to the
$L^{2}(0,1)$-norm is proportional to $(1/\nu)$ for all sufficiently small
$\nu>0$. In particular, there exists $\varrho>0$ such that for each $\nu>0$,
this $\nu$-covering number does not exceed $2^{(\varrho/2)(1/\nu)}$. Fix
$\nu>0$. Let $y_{1},\ldots,y_{K}$ be a $(\nu/\sqrt{m})$-cover of
$\overline{\mathbb{B}}_{W_{1}^{2}(0,1)}(0;1)$ with respect to the
$L^{2}(0,1)$-norm. We have $K\leq 2^{(\varrho/2)(\sqrt{m}/\nu)}$. We consider
$\mathcal{W}_{K}\coloneqq\\{w\in L^{2}(0,1;\mathbb{R}^{m})\ \colon
w_{i}\in\\{y_{1},\ldots,y_{K}\\},\quad i=1,\ldots,m\\}$. Fix
$v\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1,\mathbb{R}^{m})}(0;1)$. Hence
$v_{i}\in\overline{\mathbb{B}}_{W_{1}^{2}(0,1)}(0;1)$, $i=1,\ldots,m$. We can
choose $w\in\mathcal{W}_{K}$ such that for each $i\in\\{1,\ldots,m\\}$,
$\|v_{i}-w_{i}\|_{L^{2}(0,1)}\leq\nu/\sqrt{m}$. Hence

$\displaystyle\|v-w\|_{L^{2}(0,1;\mathbb{R}^{m})}^{2}=\sum_{i=1}^{m}\|v_{i}-w_{i}\|_{L^{2}(0,1)}^{2}\leq\nu^{2}.$

The number of elements of $\mathcal{W}_{K}$ does not exceed $K^{m}\leq
2^{(\varrho/2)m(\sqrt{m}/\nu)}$.

Let $Y_{0}=\overline{\mathbb{B}}_{W_{1}^{2}(0,1,\mathbb{R}^{m})}(0;1)$ and let
$Y=L^{2}(0,1;\mathbb{R}^{m})$. From the proof of Lemma 8.2-2 (c) in [9], we
deduce $\mathcal{N}(Y_{0};\nu)\leq\mathcal{N}(Y,Y_{0};\nu/2)$ for all $\nu>0$.
Combining the pieces ensures the assertion. ∎

##  Appendix C Uniform expectation bounds in Hilbert spaces

We establish a uniform expectation bound for sample averages of independent
Hilbert space-valued random variables using a chaining argument (see, e.g.,
[6, 27]) and an exponential moment bound [20]. The uniform expectation bound
allows us to establish convergence rates for SAA optimal values and SAA
critical points, as demonstrated in the main text.

Throughout the section, let $(\Theta,\mathcal{A},\mathbb{P})$ be a complete
probability space, and let $\xi$ be a $\Xi$-valued random element, where $\Xi$
is a complete separable metric space.

######  Theorem C.1.

Let $H$ be a real, separable Hilbert space, let $C\subset H$ be a nonempty,
closed set with diameter $D>0$, and let $G\colon C\times\Xi\to H$ be a
Carathéodory mapping with $\mathbb{E}[\|G(x,\xi)\|_{H}]<\infty$ and
$\mathbb{E}[G(x,\xi)]=0$ for each $x\in C$. Suppose there exists a constant
$M>0$ such that

$\displaystyle\mathbb{E}[\cosh(\lambda\|G(x,\xi)-G(y,\xi)\|_{H})]\leq\exp(M^{2}\lambda^{2}\|x-y\|_{H}^{2}/2)\;\text{for
all}\;x,y\in C,\;\lambda\geq 0.$ (21)

Let $\xi^{1},\xi^{2},\ldots$ be independent $\Xi$-valued random elements
defined on $(\Theta,\mathcal{A},\mathbb{P})$, each having the same
distribution as $\xi$.

We define $\widehat{G}_{N}(x)\coloneqq(1/N)\sum_{i=1}^{N}G(x,\xi^{i})$. Then
for each $x_{0}\in C$ and $N\in\mathbb{N}$,

$\displaystyle\mathbb{E}[\sup_{x\in
C}\,\|\widehat{G}_{N}(x)\|_{H}]\leq\mathbb{E}[\|\widehat{G}_{N}(x_{0})\|_{H}]+\frac{4\sqrt{3}M}{\sqrt{N}}\int_{0}^{D/2}\sqrt{\ln(2\mathcal{N}(C,\varepsilon))}\,\mathrm{d}\varepsilon.$

Under the hypotheses of Theorem C.1, we have (cf. [29, p. 79])

$\displaystyle\mathbb{E}[\|\widehat{G}_{N}(x_{0})\|_{H}]\leq\sqrt{\frac{\mathbb{E}[\|G(x_{0},\xi)\|_{H}^{2}]}{N}}.$

We leverage this bound in the main text.

Before we prepare our proof of Theorem C.1, we comment on the condition (21).

######  Remark C.2.

  1. (a)

For a mean-zero random variable $Z\colon\Xi\to\mathbb{R}$, the condition,
$\mathbb{E}[\cosh(\lambda|Z|)]\leq\exp(\lambda^{2}\tau^{2}/2)$ for all
$\lambda\in\mathbb{R}$ and some $\tau>0$, is equivalent to $Z$ being sub-
Gaussian (cf. [18, Lem. B.2]). This provides one motivation for using the
condition (21) to model sub-Gaussian behavior.

  2. (b)

If $G(\cdot,\xi)$ is Lipschitz continuous with Lipschitz constant $M$ for each
$\xi\in\Xi$, then (21) holds true. This is a consequence of the inequality
$\cosh(x)\leq\exp(x^{2}/2)$ valid for all $x\in\mathbb{R}$.

We prepare our proof of Theorem C.1. The following lemma provides a basic
upper bound on the expected pointwise maximum of sub-Gaussian random
variables.

######  Lemma C.3 (see [18, Lem. B.5]).

Let $\sigma>0$. If $Z_{k}:\Theta\to\mathbb{R}$ are random variables with
$\mathbb{E}[\cosh(\lambda|Z_{k}|)]\leq\exp(\lambda^{2}\sigma^{2}/2)$ for all
$\lambda\in\mathbb{R}$ and $k=1,2,\ldots,K$, then $\mathbb{E}[\max_{1\leq
k\leq K}|Z_{k}|]\leq\sigma\sqrt{2\ln(2K)}$.

The following proposition is inspired by entropy integral bounds for sub-
Gaussian processes (see, e.g., [5, 6]).

######  Proposition C.4.

Let $(Y_{1},d_{Y_{1}})$ be a separable metric space, let $C\subset Y_{1}$ be a
nonempty, closed set with diameter $D_{1}>0$, let $Y_{2}$ be a separable
Banach space, and let $G\colon C\times\Xi\to Y_{2}$ be a Carathéodory mapping.
Suppose that

$\displaystyle\mathbb{E}[\cosh(\lambda\|G(x,\xi)-G(y,\xi)\|_{Y_{2}})]\leq\exp(\lambda^{2}d_{Y_{1}}(x,y)^{2}/2)\;\text{for
all}\;x,y\in C,\;\lambda\geq 0.$ (22)

Then for all points $x_{0}\in C$, we have

$\displaystyle\mathbb{E}[\sup_{x\in
C}\,\|G(x,\xi)\|_{Y_{2}}]\leq\mathbb{E}[\|G(x_{0},\xi)\|_{Y_{2}}]+4\sqrt{2}\int_{0}^{D_{1}/2}\sqrt{\ln(2\mathcal{N}(C,\varepsilon))}\,\mathrm{d}\varepsilon.$
(23)

###### Proof.

We adapt the proof of Theorem 3.1 in [5, p. 95]. If the integral in (23) is
infinite, then there is nothing to show. Next, let the integral in (23) be
finite. Hence $C$ is totally bounded and $D$ is finite. We define
$\varepsilon_{0}\coloneqq D_{1}$ and
$\varepsilon_{k}\coloneqq(1/2)^{k}\varepsilon_{0}$. For each $k\in\mathbb{N}$,
let $C_{k}$ be a $\varepsilon_{k}$-net in the set $C$ with
$\mathcal{N}(C,\varepsilon_{k})$ elements. We define
$C_{\infty}\coloneqq\cup_{k=0}^{\infty}C_{k}$. The set $C_{0}$ contains only
one point. We choose this point to be $x_{0}$.

Fix $k\in\mathbb{N}\cup\\{0\\}$. Following [5, p. 94], we define a mapping
$\pi_{k}\colon C\to C_{k}$ such that $\pi_{k}(x)=x$ if $x\in C_{k}$ and
$\pi_{k}(x)$ equals a point in $C_{k}$ closest to $x$ if $x\not\in C_{k}$.

Let $x\in C_{\infty}$. Then there exists $K(x)\in\mathbb{N}\cup\\{0\\}$ such
that $x\in C_{K(x)}$. We define $x_{K(x)}\coloneqq x$, and
$x_{k-1}\coloneqq\pi_{k-1}(x_{k})$, for $k=1,\ldots,K(x)$. We have

$\displaystyle
G(x,\xi)=G(x_{0},\xi)+\sum_{k=1}^{K(x)}(G(x_{k},\xi)-G(x_{k-1},\xi)).$

Hence

$\displaystyle\|G(x,\xi)\|_{Y_{2}}\leq\|G(x_{0},\xi)\|_{Y_{2}}+\sum_{k=1}^{K(x)}\max_{x\in
C_{k}}\|G(x,\xi)-G(\pi_{k-1}(x),\xi)\|_{Y_{2}}.$

Combined with (22), $d_{Y_{1}}(x,\pi_{k-1}(x))\leq\varepsilon_{k-1}$ for all
$x\in C$, and Lemma C.3, we find that for $k=1,2,\ldots$,

$\displaystyle\mathbb{E}[\max_{x\in
C_{k}}\,\|G(x,\xi)-G(\pi_{k-1}(x),\xi)\|_{Y_{2}}]\leq\varepsilon_{k-1}\sqrt{2\ln(2\mathcal{N}\big{(}C,\varepsilon_{k}\big{)})}.$

We obtain

$\displaystyle\mathbb{E}[\sup_{x\in C_{\infty}}\,\|G(x,\xi)\|_{Y_{2}}]$
$\displaystyle\leq\|G(x_{0},\xi)\|_{Y_{2}}+\sum_{k=1}^{\infty}\,\varepsilon_{k-1}\sqrt{2\ln(2\mathcal{N}\big{(}C,\varepsilon_{k}\big{)})}.$

Since the set $C_{\infty}$ is dense in $C$, $G$ is a Carathéodory mapping, and

$\displaystyle\varepsilon_{k-1}\sqrt{2\ln(2\mathcal{N}\big{(}C,\varepsilon_{k}\big{)})}\leq
4\int_{\varepsilon_{k+1}}^{\varepsilon_{k}}\sqrt{2\ln(2\mathcal{N}\big{(}C,\varepsilon\big{)})}\mathrm{d}\varepsilon,$

we obtain the assertion. ∎

The next lemma provides an exponential moment inequality. It is a direct
consequence of Theorem 3 in [20] and Lemma 1 in [14], and has been established
in [18].

######  Lemma C.5 (see [18, Lem. B.4]).

Let $\tau>0$, let $H$ be a real, separable Hilbert space, and let
$Z_{i}:\Theta\to H$ be independent, mean-zero random vectors such that
$\mathbb{E}[\cosh(\lambda\|Z_{i}\|_{H})]\leq\exp(\lambda^{2}\tau^{2}/2)$ for
all $\lambda\geq 0$ and $i\in\\{1,2,\ldots,N\\}$. Then for each $\lambda\geq
0$,

$\mathbb{E}[\cosh(\lambda\|Z_{1}+\cdots+Z_{N}\|_{H})]\leq\exp(3\lambda^{2}\tau^{2}N/4).$

Now we are ready to establish Theorem C.1.

###### Proof of Theorem C.1.

Defining the metric
$d_{Y_{1}}(x,y)\coloneqq\sqrt{3/2}(M/\sqrt{N})\|x-y\|_{H}$, Propositions C.4
and C.5 imply the assertion. ∎

The following example demonstrates an upper bound on the integral in (23)
under a bound on the covering numbers. We use it in the main text to establish
mean convergence rates for SAA optimal values and criticality measures.

######  Example C.6.

Let $c_{1}>0$ be a constant, and let $C\subset Y_{1}$ be a subset of a metric
space $(Y_{1},d_{Y_{1}})$ with diameter $D>0$. If
$\mathcal{N}(C,\varepsilon)\leq 2^{c_{1}/\varepsilon}$ for all
$\varepsilon>0$, then

$\displaystyle\int_{0}^{D/2}\sqrt{\ln(2\mathcal{N}(C,\varepsilon))}\,\mathrm{d}\varepsilon\href
https://wolfr.am/1nYJyzwqF\sqrt{D/2}(\sqrt{c_{1}+D/2}+\sqrt{c_{1}}).$

## References

  * [1] L. Ambrosio and G. Dal Maso.  A general chain rule for distributional derivatives.  Proc. Amer. Math. Soc., 108(3):691–702, 1990.  doi:10.2307/2047789.
  * [2] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl.  CasADi: a software framework for nonlinear optimization and optimal control.  Math. Program. Comput., 11(1):1–36, 2019.  doi:10.1007/s12532-018-0139-4.
  * [3] H. H. Bauschke and P. L. Combettes.  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  CMS Books Math. Springer, New York, 2011.  doi:10.1007/978-1-4419-9467-7.
  * [4] M. Š. Birman and M. Z. Solomjak.  Quantitative Analysis in Sobolev Imbedding Theorems and Applications to Spectral Theory.  Amer. Math. Soc. Transl. Ser. 2, 114. AMS, Providence, RI, 1980.  Translated from the Russian by F. A. Cezus.  doi:10.1090/trans2/114.
  * [5] V. V. Buldygin and Yu. V. Kozachenko.  Metric Characterization of Random Variables and Random Processes.  Transl. Math. Monogr. 188. AMS, Providence, RI, 2000.  Translated by V. Zaiats.  doi:10.1090/mmono/188.
  * [6] E. Giné and R. Nickl.  Mathematical Foundations of Infinite-Dimensional Statistical Models.  Camb. Ser. Stat. Probab. Math. 40. Cambridge University Press, New York, 2016.  doi:10.1017/CBO9781107337862.
  * [7] T. Hytönen, J. van Neerven, M. Veraar, and L. Weis.  Analysis in Banach Spaces: Martingales and Littlewood-Paley Theory.  Ergeb. Math. Grenzgeb. (3) 63. Springer, Cham, 2016.  doi:10.1007/978-3-319-48520-1.
  * [8] A. J. Kleywegt, A. Shapiro, and T. Homem-de Mello.  The sample average approximation method for stochastic discrete optimization.  SIAM J. Optim., 12(2):479–502, 2002.  doi:10.1137/S1052623499363220.
  * [9] E. Kreyszig.  Introductory Functional Analysis with Applications.  John Wiley & Sons, New York, NY, 1978.
  * [10] S. Lenhart and J. T. Workman.  Optimal control applied to biological models.  Chapman & Hall/CRC Mathematical and Computational Biology Series. Chapman & Hall/CRC, Boca Raton, FL, 2007.  doi:10.1201/9781420011418.
  * [11] G. Leoni.  A first course in Sobolev spaces, volume 181 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2nd edition, 2017.  doi:10.1090/gsm/181.
  * [12] T. Lew, R. Bonalli, and M. Pavone.  Sample average approximation for stochastic programming with equality constraints, 2022.  doi:10.48550/ARXIV.2206.09963.
  * [13] R. Miller Neilan and S. Lenhart.  An introduction to optimal control with an application in disease modeling.  In Modeling paradigms and analysis of disease transmission models, volume 75 of DIMACS Ser. Discrete Math. Theoret. Comput. Sci., pages 67–81. Amer. Math. Soc., Providence, RI, 2010.  doi:10.1090/dimacs/075/03.
  * [14] J. Milz.  Sample average approximations of strongly convex stochastic programs in Hilbert spaces.  Optim. Lett., 17:471–492, 2023.  doi:10.1007/s11590-022-01888-4.
  * [15] J. Milz.  EnsembleControl, June 2024.  doi:10.5281/zenodo.11669862.
  * [16] J. Milz.  Supplementary code for the manuscript: Convergence rates for ensemble-based solutions to optimal control of uncertain dynamical systems, July 2024.  doi:10.5281/zenodo.12740932.
  * [17] J. Milz and T. M. Surowiec.  Asymptotic consistency for nonconvex risk-averse stochastic optimization with infinite dimensional decision spaces.  Math. Oper. Res., 2023.  doi:10.1287/moor.2022.0200.
  * [18] J. Milz and M. Ulbrich.  Sample size estimates for risk-neutral semilinear PDE-constrained optimization.  SIAM J. Optim., 34(1):844–869, 2024.  doi:10.1137/22M1512636.
  * [19] C. Phelps, J. Royset, and Q. Gong.  Optimal control of uncertain systems using sample average approximations.  SIAM J. Control Optim., 54(1):1–29, 2016.  doi:10.1137/140983161.
  * [20] I. F. Pinelis and A. I. Sakhanenko.  Remarks on inequalities for large deviation probabilities.  Theory Probab. Appl., 30(1):143–148, 1986.  doi:10.1137/1130013.
  * [21] E. Polak.  Optimization: Algorithms and Consistent Approximations.  Appl. Math. Sci. 124. Springer, New York, 1997.  doi:10.1007/978-1-4612-0663-7.
  * [22] W. Römisch and T. M. Surowiec.  Asymptotic properties of Monte Carlo methods in elliptic PDE-constrained optimization under uncertainty.  preprint, https://arxiv.org/abs/2106.06347, 2021.
  * [23] P. T. Roy, A. B. Owen, M. Balandat, and M. Haberland.  Quasi-Monte Carlo methods in Python.  Journal of Open Source Software, 8(84):5309, 2023.  doi:10.21105/joss.05309.
  * [24] D. Ruppen, C. Benthack, and D. Bonvin.  Optimization of batch reactor operation under parametric uncertainty — computational aspects.  Journal of Process Control, 5(4):235–240, 1995.  doi:10.1016/0959-1524(95)00015-I.
  * [25] A. Scagliotti.  Optimal control of ensembles of dynamical systems.  ESAIM Control Optim. Calc. Var., 29:Paper No. 22, 39, 2023.  doi:10.1051/cocv/2023011.
  * [26] B. Srinivasan, D. Bonvin, E. Visser, and S. Palanki.  Dynamic optimization of batch processes: II. role of measurements in handling uncertainty.  Computers & Chemical Engineering, 27(1):27–44, 2003.  doi:10.1016/S0098-1354(02)00117-5.
  * [27] M. Talagrand.  Upper and lower bounds for stochastic processes—decomposition theorems, volume 60 of Ergebnisse der Mathematik und ihrer Grenzgebiete. 3. Folge. A Series of Modern Surveys in Mathematics.  Springer, Cham, 2nd edition, 2021.  doi:10.1007/978-3-030-82595-9.
  * [28] A. Wächter and L. T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Math. Program., 106(1):25–57, 2006.  doi:10.1007/s10107-004-0559-y.
  * [29] V. Yurinsky.  Sums and Gaussian Vectors.  Lecture Notes in Math. 1617. Springer, Berlin, 1995.  doi:10.1007/BFb0092599.


# Biomedical Entity Linking with Triple-aware Pre-Training

Xi Yan
HITEC e.V.
Universität Hamburg
<EMAIL_ADDRESS>
Cedric Möller
Universität Hamburg
<EMAIL_ADDRESS>
Ricardo Usbeck
Universität Hamburg
<EMAIL_ADDRESS>

###### Abstract

Linking biomedical entities is an essential aspect in biomedical natural
language processing tasks, such as text mining and question answering.
However, a difficulty of linking the biomedical entities using current large
language models (LLM) trained on a general corpus is that biomedical entities
are scarcely distributed in texts and therefore have been rarely seen during
training by the LLM. At the same time, those LLMs are not aware of high level
semantic connection between different biomedical entities, which are useful in
identifying similar concepts in different textual contexts. To cope with
aforementioned problems, some recent works Yuan et al. (2022b); Singh et al.
(2022) focused on injecting knowledge graph information into LLMs. However,
former methods either ignore the relational knowledge of the entities or lead
to catastrophic forgetting. Therefore, we propose a novel framework to pre-
train the powerful generative LLM by a corpus synthesized from a KG. In the
evaluations we are unable to confirm the benefit of including synonym,
description or relational information.

##  1 Introduction

Biomedical entity linking (EL) is a critical process in biomedical text mining
that seeks to identify and associate relevant biological and medical entities
mentioned in unstructured text with their corresponding identifiers in
knowledge bases. Accurate recognition and linking of these entities are
pivotal in promoting biomedical research, drug discovery, and personalized
medicine Chandak et al. (2023). Although substantial progress has been made in
recent years, there is an ongoing need for refining methods and techniques
employed for entity linking in the biomedical domain.

In this report, we present a novel approach that integrates linearized triples
into the biomedical entity linking process while reevaluating the inclusion of
synonym information. Our proposed method linearizes triples and considers them
during the pre-training step. In past studies, synonym information, which
involves using alternative names or terminologies for the same biomedical
entity, has been suggested as a means to enhance entity linking when used
during pre-training Yuan et al. (2022b); Xu et al. (2023). Our study aims to
build upon this existing knowledge by integrating both strategies and
assessing their impact on performance.

Despite the reported benefits of synonym information in prior studies, our
analysis of this approach, combined with the introduction of linearized
triples Li et al. (2021), yielded different results. We find that
incorporating linearized triples only lead to minimal improvements in our
entity linking model’s performance. Moreover, we are unable to confirm the
purported advantages of including synonym information in our experiments,
which stands in contrast to the findings of previous literature.

We highlight the limitations of our study and suggest possible avenues for
future research to further advance biomedical entity linking techniques by
building on our work with linearized triples and reevaluating synonym
information.

##  2 Related work

Entity Linking has a long history of research. Recent methods can be
categorized into two types. First, discriminative methods that are based on
the bi-encoder / cross-encoder pairing Wu et al. (2020a); Logeswaran et al.
(2019); Ayoola et al. (2022). Both encoders are commonly BERT-like models. The
bi-encoder encodes the description of each entity and matches it to the text
by using an approximate nearest neighbor search. This is important as the next
step, the cross-encoding, is expensive. Here, those neighbors are reranked by
applying a cross-encoder to the concatenation of both, the input text and the
entity description. The highest-ranked entity is then the final linked one. In
the biomedical domain, the works by Angell et al. (2021), Varma et al. (2021),
Agarwal et al. (2021) and Bhowmik et al. (2021) fall into this category.

Another type of entity linker is based on generative models Cao et al. (2021,
2022); Xu et al. (2023). Here, instead of using some external description of
an entity, the whole model memorizes the KG during training. The linked entity
is then directly generated by the model. Such methods skip the problem of
mining negatives which are crucial for a good performance of bi-encoder-based
methods. Only the work by Yuan et al. (2022b) is based on such methods in the
biomedical domain. As generative models lack the ability to incorporate
external information, they alleviate this problem by introducing a pre-
training stage where syntactical information from a knowledge graph is
learned. This is especially important in the biomedical domain as entities
often own a large variety of synonyms. We build upon their work by extending
the pre-training regime to the inclusion of triple information.

##  3 Method

###  3.1 Task definition

Given are a text $t$, a set of marked mentions $M_{t}$ in the text and a KG
$\mathcal{G}=(\mathcal{E},\mathcal{R},E)$. The KG consists of a set of
entities $\mathcal{E}$, a set of relations $\mathcal{R}$ and a set of edges
composed of head entity, relation and tail entity
$E\subseteq(\mathcal{E}\times\mathcal{R}\times\mathcal{E})$. The task is to
identify the subset of entities $E_{t}\subseteq\mathcal{E}$ which the mentions
$M_{t}$ are referring to.

###  3.2 Model

In the vein of the work by Cao et al. (2021), we model the problem as a
sequence-to-sequence generation task. The input to the generative model is
text and the output are the generated entity identifiers in the corresponding
KGs. Similar to other works Cao et al. (2021, 2022); Yuan et al. (2022b), we
consider the definition of the concepts in the corresponding KGs as the unique
textual representation of each concept.

###  3.3 Pre-training

We linearize the information from synonym and triples in the pre-training
stage. They are linearized into a synthesized corpora before feeding into the
BART. We have tested 2 different settings for converting the triples, namely
line-by-line and all-in-one.

In terms of the synonym information, we follow the setting by Yuan et al.
(2022b). We first extract the description of the entity and convert it to a
text of the following form:

⬇

[BOS][ST]$s^{a}_{e}$[ET] is defined as $c_{e}$[EOS]

Here, $s^{a}_{e}$ stands for the synonym $a$ and $c_{e}$ for the description
of entity $e$. This would be the input to the encoder of the generative model.

As an output the model has to generate:

⬇

[BOS] $s_{e}^{a}$ is $s_{e}^{b}$ [EOS]

This lets the model learn the connection between the different synonyms of the
same entity.

Based on that, we introduce an additional pre-training step to incorporate
more semantic information by utilising triple information from the underlying
knowledge graph. A triple is of the form $<e,r,e^{\prime}>$ which describes
that a relationship $r$ holds between entity $e$ and $e^{\prime}$. The input
is here the same as for the synonym information. The output is of the form:

⬇

[BOS] $s_{e}^{a}$ $l_{r}$ $s^{b}_{e^{\prime}}$ [EOS]

$l_{r}$ is here the label of relation $r$. We denote this line-by-line

Furthermore, we experimented with an all-in-one pre-training approach of the
form:

⬇

[BOS] $s^{a}_{e}$ $l_{r_{1}}$ $s^{b_{1}}_{e_{1}}$ $\dots$ $l_{r_{n}}$
$s^{b_{n}}_{e_{n}}$ [EOS]

See Figure 1 for an overview of the pre-training.

Figure 1: An overall workflow of our framework. We adopt different
textualization formats for synonym information and triples. Both are included
in the pre-training stage.

###  3.4 Fine-tuning

During fine-tuning, the model is trained for the actual entity linking task.
The input to the generative model is the unlabelled biomedical text. To
generate the linked entities, each mention is included in a template as
follows:

⬇

[BOS] $m_{i}$ is $s^{a}_{e}$ [EOS]

The model then generates the entity identifier after the token "is". Similar
to the work by Yuan et al. (2022b), we choose the synonym which is
syntactically close to the corresponding mention in the text as the target
entity identifier during fine-tuning.

The generated entity identifier is mapped back to the concrete entity in the
final step via a lookup table. During inference, we restrict the possible
output space by limiting it to the available entity names and synonyms.

See Figure 2 for an overview of the pre-training.

Figure 2:  An overview of the fine-tuning stage

##  4 Evaluation

###  4.1 Pre-training Strategy

We use a synthesized corpus composed of triples, synonyms and descriptions
from UMLS. More specifically, we decide to use a subset of UMLS, i.e., st21pv
Mohan and Li (2019). It is a well-connected KG with information about concept
definitions and synonyms. Specifically, 160K out of 2.37M concepts have
definitions, 1.11M concepts have several synonyms and 68K concepts are
connected to on average 8 triples. During the pre-training step, we construct
samples by iterating through each concept’s synonyms and triples. Each concept
is densely connected and the distribution of the number of triples a concept
is connected to is skewed. For instance, some "popular" concepts are connected
to over 1000 triples, while some are connected to only 1 triple. To avoid the
class imbalance, we sample the included triples based on the relation
frequencies.

To train the model with KG information, we linearize triples. Linearization
refers to a special type of technique on converting graph to text, i.e.,
converting triples to one/more sentences which serve as input of the LLM.

We sample the included triples based on the relation frequencies. First, we
gather the occurrence frequency of all relations in the KB by counting the
number of triples this relation is connected to. We define the probability
($P_{r}$) of a relation $r$ to be negatively related to the frequency. Then,
for each concept in the KG, we collect its connected triples and segment the
triples into different groups based on their relation $r$.

Both settings are trained under the same experiment setting with a batch size
of 128. We save the best model within 12 training epochs.

####  4.1.1 Fine-tuning

The model is fine-tuned on two established datasets, namely BC5CDR (Li et al.
(2016)) and NCBI (Dogan et al. (2014)). Those entity linking datasets are
constructed on subsets of UMLS, making them perfect choices to test our
model’s performance on. Among the datasets, NCBI and BC5CDR are generated by
annotating PubMed papers. On the other hand, NCBI and BC5CDR are annotated
against Medical Subject Headings (MeSH) - a terminology knowledge graph for
indexing and cataloging of biomedical information.

The statistics of the four datasets are exhibited in Table 1 below. As we can
see, NCBI and BC5CDR (annotated on academic text) are smaller in size. Also
NCBI and BC5CDR are dense in terms of the target entities they contain (14,967
and 268,162).

Nums | NCBI | BC5CDR
---|---|---
Train | 5,784 | 9,285
Dev | 787 | 9,515
Test | 960 | 9,654
Entities | 14,967 | 268,162
Table 1: Numbers of the samples in the training, development and test set

BART-large Liu et al. (2023) is chosen as the generative model as it has been
an established benchmark model for such tasks.

###  4.2 Results

We assess the performance of four distinct models in the entity linking task,
including two of our own models, each pre-trained via either a line-by-line or
all-in-one strategy, a synonym pre-trained model from Yuan et al. (2022b)
(denoted Syn-Only), and a basic BART model. We also include the recent papers
which pretrains BART on biomedical domain Yuan et al. (2022a) before finetuned
on biomedical entity linking datasets and ResCNN Lai et al. (2021) which
achieves state-of-the-art results on various biomedical EL datasets. Each
model undergoes fine-tuning specific to the entity linking task. Recall@1 for
each model are presented in Tables LABEL:table:result1 and 2. We limit
ourselves to Recall@1 to follow the common practice when measuring entity
linking performance without named entity recognition. The best-performing
metrics are emphasized in bold.

| BC5CDR | NCBI
---|---|---
Syn-Only | 93.3% | 91.9%
Syn-Only | 92.68% | 89.45%
All-in-one | 92.86% | 88.43%
Line-by-line | 92.66% | 90.00%
BART | 92.58% | 89.06%
BioBART-Large | 93.01% | 89.27%
BioBART-Base |  93.26% | 89.40%
ResCNN | 91.7 % |  92.4%
Table 2: Recall@1 on BC5CDR and NCBI, which are PubMed articles annotated
against MESH.

###  4.3 Analysis

Based on the table, our triple injection framework exceeds the BART baseline
on the 2 Wbenchmarks datasets. On BC5CDR and NCBI, the gain compared to BART
is around 0.2% and 0.5%.

Does triple injection enhance model’s capacity to link to the correct entity?
The answer is yes, since over 4 datasets, the All-in-one or Line-by-line
variants outperform the variant that was not trained on the linearized corpora
for around 1% (Recall@1).

##  5 Conclusion

Our study seek to improve biomedical entity linking through the integration of
linearized triples and synonym information. However, contrary to expectation,
the incorporation of these elements leads to only minimal improvements in our
EL model performance.

In conclusion, our study underscores the complexities of biomedical EL and
prompts the need for more sophisticated approaches to improve its accuracy. A
possible future extension of this work could be to explore more sophisticated
methods to instruct the LLMs to learn external knowledge, such that the
knowledge is injected in an efficient way which benefits the models in
downstream tasks. For instance, by incorporating the KG information not just
in a linearized manner but by exploiting the graph-structure with Graph Neural
Networks Wu et al. (2020b), mutliple methods could be further developed.

This work has been partially supported by grants for the DFG project
NFDI4DataScience project (DFG project no. 460234259) and by the Federal
Ministry for Economics and Climate Action in the project CoyPu (project number
01MK21007G).

## References

  * Agarwal et al. (2021) Dhruv Agarwal, Rico Angell, Nicholas Monath, and Andrew McCallum. 2021.  Entity linking and discovery via arborescence-based supervised clustering.  _CoRR_ , abs/2109.01242.
  * Angell et al. (2021) Rico Angell, Nicholas Monath, Sunil Mohan, Nishant Yadav, and Andrew McCallum. 2021\.  Clustering-based inference for biomedical entity linking.  In _Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL-HLT 2021, Online, June 6-11, 2021_ , pages 2598–2608. Association for Computational Linguistics.
  * Ayoola et al. (2022) Tom Ayoola, Shubhi Tyagi, Joseph Fisher, Christos Christodoulopoulos, and Andrea Pierleoni. 2022.  Refined: An efficient zero-shot-capable approach to end-to-end entity linking.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies: Industry Track, NAACL 2022, Hybrid: Seattle, Washington, USA + Online, July 10-15, 2022_ , pages 209–220. Association for Computational Linguistics.
  * Bhowmik et al. (2021) Rajarshi Bhowmik, Karl Stratos, and Gerard de Melo. 2021.  Fast and effective biomedical entity linking using a dual encoder.  In _Proceedings of the 12th International Workshop on Health Text Mining and Information Analysis, LOUHI@EACL, Online, April 19, 2021_ , pages 28–37. Association for Computational Linguistics.
  * Cao et al. (2021) Nicola De Cao, Gautier Izacard, Sebastian Riedel, and Fabio Petroni. 2021.  Autoregressive entity retrieval.  In _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net.
  * Cao et al. (2022) Nicola De Cao, Ledell Wu, Kashyap Popat, Mikel Artetxe, Naman Goyal, Mikhail Plekhanov, Luke Zettlemoyer, Nicola Cancedda, Sebastian Riedel, and Fabio Petroni. 2022.  Multilingual autoregressive entity linking.  _Trans. Assoc. Comput. Linguistics_ , 10:274–290.
  * Chandak et al. (2023) Payal Chandak, Kexin Huang, and Marinka Zitnik. 2023.  Building a knowledge graph to enable precision medicine.  _Scientific Data_ , 10(1):67.
  * Dogan et al. (2014) Rezarta Islamaj Dogan, Robert Leaman, and Zhiyong Lu. 2014.  NCBI disease corpus: A resource for disease name recognition and concept normalization.  _J. Biomed. Informatics_ , 47:1–10.
  * Lai et al. (2021) Tuan Lai, Heng Ji, and ChengXiang Zhai. 2021.  BERT might be overkill: A tiny but effective biomedical entity linker based on residual convolutional neural networks.  In _Findings of the Association for Computational Linguistics: EMNLP 2021_ , pages 1631–1639, Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Li et al. (2016) Jiao Li, Yueping Sun, Robin J. Johnson, Daniela Sciaky, Chih-Hsuan Wei, Robert Leaman, Allan Peter Davis, Carolyn J. Mattingly, Thomas C. Wiegers, and Zhiyong Lu. 2016.  Biocreative V CDR task corpus: a resource for chemical disease relation extraction.  _Database J. Biol. Databases Curation_ , 2016.
  * Li et al. (2021) Junyi Li, Tianyi Tang, Wayne Xin Zhao, Zhicheng Wei, Nicholas Jing Yuan, and Ji-Rong Wen. 2021.  Few-shot knowledge graph-to-text generation with pretrained language models.  In _Findings of the Association for Computational Linguistics: ACL-IJCNLP 2021_ , pages 1558–1568, Online. Association for Computational Linguistics.
  * Liu et al. (2023) Pengfei Liu, Weizhe Yuan, Jinlan Fu, Zhengbao Jiang, Hiroaki Hayashi, and Graham Neubig. 2023.  Pre-train, prompt, and predict: A systematic survey of prompting methods in natural language processing.  _ACM Computing Surveys_ , 55(9):1–35.
  * Logeswaran et al. (2019) Lajanugen Logeswaran, Ming-Wei Chang, Kenton Lee, Kristina Toutanova, Jacob Devlin, and Honglak Lee. 2019.  Zero-shot entity linking by reading entity descriptions.  In _Proceedings of the 57th Conference of the Association for Computational Linguistics, ACL 2019, Florence, Italy, July 28- August 2, 2019, Volume 1: Long Papers_ , pages 3449–3460. Association for Computational Linguistics.
  * Mohan and Li (2019) Sunil Mohan and Donghui Li. 2019.  Medmentions: A large biomedical corpus annotated with UMLS concepts.  In _1st Conference on Automated Knowledge Base Construction, AKBC 2019, Amherst, MA, USA, May 20-22, 2019_.
  * Singh et al. (2022) Ayush Singh, Saranya Krishnamoorthy, and John Ortega. 2022.  Neighbert: Medical entity linking using relation induced dense retrieval.
  * Varma et al. (2021) Maya Varma, Laurel J. Orr, Sen Wu, Megan Leszczynski, Xiao Ling, and Christopher Ré. 2021.  Cross-domain data integration for named entity disambiguation in biomedical text.  In _Findings of the Association for Computational Linguistics: EMNLP 2021, Virtual Event / Punta Cana, Dominican Republic, 16-20 November, 2021_ , pages 4566–4575. Association for Computational Linguistics.
  * Wu et al. (2020a) Ledell Wu, Fabio Petroni, Martin Josifoski, Sebastian Riedel, and Luke Zettlemoyer. 2020a.  Scalable zero-shot entity linking with dense entity retrieval.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing, EMNLP 2020, Online, November 16-20, 2020_ , pages 6397–6407. Association for Computational Linguistics.
  * Wu et al. (2020b) Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip. 2020b.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1):4–24.
  * Xu et al. (2023) Zhenran Xu, Yulin Chen, and Baotian Hu. 2023.  Improving biomedical entity linking with cross-entity interaction.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 37(11):13869–13877.
  * Yuan et al. (2022a) Hongyi Yuan, Zheng Yuan, Ruyi Gan, Jiaxing Zhang, Yutao Xie, and Sheng Yu. 2022a.  BioBART: Pretraining and evaluation of a biomedical generative language model.  In _Proceedings of the 21st Workshop on Biomedical Language Processing_ , pages 97–109, Dublin, Ireland. Association for Computational Linguistics.
  * Yuan et al. (2022b) Hongyi Yuan, Zheng Yuan, and Sheng Yu. 2022b.  Generative biomedical entity linking via knowledge base-guided pre-training and synonyms-aware fine-tuning.  In _Proceedings of the 2022 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, NAACL 2022, Seattle, WA, United States, July 10-15, 2022_ , pages 4038–4048. Association for Computational Linguistics.