Convex operator-theoretic methods in stochastic control

The systematic mathematical study of stochastic systems goes back to the work of Kolmogorov [30], although related continuous-time electro-dynamic, molecular, quantum-theoretic, and thermal diffusions have been studied even earlier by Fokker [25], von Smoluchowski [45], Planck [37], and Chapman [20]. The relevance of such diffusions in the context of analyzing the ergodic behavior of dynamical systems has been observed by Has’minskiǐ [27] in 1960, who found an elegant way of using Harnack’s inequality in combination with a weak Lyapunov function in order to analyze the existence of steady states of FPK equations. And, in fact, it was only six years later that Wonham realized the great potential of Has’minskiǐ’s work in the context of ergodic control theory [46]. By now, these historical developments are well-understood and have been extended by many authors, as reviewed in great detail in [15] focusing on FPK equations and in [33] focusing on controlled diffusion processes.

This paper presents novel results about the existence and uniqueness of solutions to ergodic and infinite horizon stochastic optimal control problems. They are summarized, respectively, in Theorem 1 and Theorem 2.

In order to put the above contributions into context, this paper is structured as follows.

• •

  Section 2 reviews the relations between controlled McKean-Vlasov systems and FPKs in the presence of reflecting boundary conditions.

• •

  Section 3 explains how to use energy functionals to establish dissipativity results for optimally controlled FPKs, leading to the first main result on ergodic optimal control that is summarized in Theorem 1.

• •

  Section 4 proposes a novel approach to analyzing parabolic HJBs with Neumann boundary conditions in a suitably weighted Hilbert-Sobolev space.

• •

  Section 5 introduces Lyapunov diffusion operators in order to analyze the stabilizability of controlled diffusions via a generalized Bakry-Emery condition.

• •

  Section 6 introduces ergodic HJBs and presents the main result on infinite horizon optimal control that is summarized in Theorem 2.

• •

  Section 7 presents a finite element Galerkin projection based operator-theoretic dynamic programming method for solving parabolic and ergodic HJBs. A case study illustrates its performance for an infinite horizon stochastic optimal control problem with indefinite stage costs and non-trivial ergodic behavior.

Let $\Omega\subseteq\mathbb{R}^{\mathsf{n}}$ denote an open set, $\overline{\Omega}=\mathrm{cl}(\Omega)$ its closure and $\partial\Omega=\overline{\Omega}\setminus\Omega$ its boundary. The sets

denote—as usual—the set of functions that are, respectively, locally $L^{k}$-integrable, $L^{k}$-integrable, smooth with compact support, $k$-times continuously differentiable with continuous extensions to the boundary, and $L^{2}$-integrable with $L^{2}$-integrable weak derivatives on $\Omega$. The inner products

are defined such that $L^{2}(\Omega)$ and $H^{1}(\Omega)$ are Hilbert spaces. Here, $\nabla$ denotes the gradient operator. Additionally, if $\omega\in L_{\mathrm{loc}}^{1}(\Omega)$ with $\omega>0$ is a strictly positive weight function, the shorthands

are used. The associated inner products are then, of course, weighted with $\omega$, too. For given time intervals $[0,T]\subseteq\mathbb{R}$ and any Banach space $(X,\|\cdot\|_{X})$, we use the symbol $L^{k}(0,T;X)$ to denote the set of strongly measurable functions $\varphi:[0,T]\to X$ with

such that $L^{k}(0,T;X)$ is a Bochner space. Weak derivatives with respect to time are indicated by a “dot”; and $H^{1}(\omega)^{*}$ denotes the dual of $H^{1}(\omega)$. The space

is then again a Hilbert space with respect to its natural inner product [28, Thm. 1.31]. Finally, the natural duality pairing between a Hilbert space and its dual is denoted by $\langle\cdot,\cdot\rangle$ assuming that it is clear from the context which Hilbert space is meant.

Deterministic control-affine systems have the form

Here, $\mathbb{T}\subseteq\mathbb{R}$ denotes an open time horizon. It is finite if $T<\infty$ and infinite if $T=\infty$. The right-hand functions $f\in C^{1}(\overline{\Omega})^{n_{x}}$ and $G\in C^{1}(\overline{\Omega})^{n_{u}}$ are assumed to be defined on the closure of an open domain $\Omega\subseteq\mathbb{R}^{n_{x}}$. Moreover, $x(t)$ denotes the state and $u(t)$ the control at time $t\in\mathbb{T}$.

Let us temporarily assume that $\Omega=\mathbb{R}^{n_{x}}$. In this case, we can add a process noise to the right-hand of (2), arriving at a stochastic differential equation (SDE) of the form

where $W_{t}$, $t\geq 0$, denotes a Wiener process, $\epsilon>0$ a diffusion constant, and $x_{0}\in\mathbb{R}^{n_{x}}$ an initial state. The state of this SDE at time $t\in\mathbb{T}$ is denoted by $X(t)$. It is a random variable, whose probability density function is denoted by $\rho(t)\in L^{1}(\Omega)$. This means that for any given Borel set $S\subseteq\Omega$, the probability of the event $X(t)\in S$ is given by

Conditions under which (3) admits a martingale solution $X$ that has such an $L^{1}$-integrable probability density function $\rho$ are discussed below.

The states of many nonlinear systems have no physical interpretation outside of their natural domain $\Omega$. For instance the concentration of the reactants in a controlled chemical process are usually required to be non-negative implying that $\Omega\neq\mathbb{R}^{n_{x}}$. In such a case, the stochastic process noise model in (3) makes no sense, as the Wiener process $W_{t}$ is unbounded. Instead, at least for the case that $\Omega$ is convex, a physically more meaningful stochastic model is given by

This is a McKean-Vlasov SDE with reflection [21, 40], where the additional drift term,

acts as reflector. We have $N_{\overline{\Omega}}(x)=\{0\}$ for $x\in\Omega$, which means that (3) and (5) coincide on $\Omega$. For $x\in\partial\Omega$, however, $N_{\overline{\Omega}}(x)$ is equal to the normal cone of $\overline{\Omega}$ at $x$. Consequently, whenever (5) admits a solution $X$ that has an $L^{1}$-integrable probability density $\rho$, we expect that

which is the case under very mild assumptions that are, for instance, discussed in [21], but sufficient conditions under which this is the case are also discussed below. For the special case that $\Omega=\mathbb{R}^{n_{x}}$, (5) reduces to (3).

Let $\mathcal{A}$ denote the $\epsilon$-diffusion operator of $f$ and $\mathcal{B}$ the linear transport operator of $G$, formally defined by

where “$\mathrm{div}$” denotes the divergence and “$\Delta$” the Laplace operator. The pair $(\mathcal{A},\mathcal{B})$ is associated with an infinite-dimensional linear control system, called the controlled FPK equation. In the presence of reflecting boundary conditions, it has the form

Here, $\rho_{0}\in L^{1}(\Omega)$ denotes a given initial probability distribution on $\Omega$, while

denote, respectively, the FPK’s open domain and its spatial boundary. Moreover, $n:\partial\Omega\to\mathbb{R}^{n_{x}}$ denotes an outer normal, such that

The open-loop linear control system (10) is closely related to the closed-loop nonlinear control system (5) after identifying $\rho$ with the probability distribution of the random process $X$ and after substituting

In order to elaborate on this relation between (5), (10) and (11), a couple of definitions are needed. For this aim, the shorthand

is introduced, where $\omega\in L_{\mathrm{loc}}^{1}(\Omega)$, $\omega>0$, is a strictly positive function. Moreover, the linear operators $\mathcal{A}$ and $\mathcal{B}$ are associated with the bilinear forms

which are defined for all $\rho\in\mathcal{P}(\omega^{-1})$ and all $V\in\mathcal{P}(\omega)$.

From a control theoretic perspective, it is interesting to point out that for the purpose of analyzing the finite-dimensional closed-loop controlled nonlinear system $(f,G)$, as defined by (5), it sufficient to analyze the infinite dimensional open-loop controlled linear system $(\mathcal{A},\mathcal{B})$, as defined in (10). This is because the existence of a unique weak solution $\rho\in\mathcal{P}(\omega^{-1})$ of (10) for $\omega\in L^{\infty}(\Omega)$ and an initial probability density function $0\leq\rho_{0}\in L^{2}(\omega^{-1})\cap L^{1}(\Omega)$ with $\int_{\Omega}\rho_{0}\mathrm{d}x=1$ of the initial random variable $X(0)$ is sufficient to ensure that (5) has a martingale solution $X$ with probability density function $\rho$ that satisfies (6). This is a standard result from the theory of diffusion processes [42], which has recently been extended for general FPKs with reflecting boundary conditions [7]; see also [21]. Of course, here we have to substitute (11), but this substitution is justified, too. Namely, it follows from the minimum principle for parabolic FPKs [15] that if $\rho\in\mathcal{P}(\omega^{-1})$, is a weak solution of (10) for an initial probability density $\rho_{0}\geq 0$, then

recalling that $\Omega_{\mathbb{T}}$ is an open set, $0\notin\mathbb{T}$, and $\epsilon>0$. Thus, the closed-loop control law, $\mu=\nu/\rho$, can be recovered from the infinite-dimensional open-loop input $\nu$ and the strictly positive function $\rho$.

Let $\mathbb{U}\subseteq\mathbb{R}^{n_{u}}$ denote a control constraint set, let

denote the associated set of feasible control laws, and let $\ell:\Omega\times\mathbb{U}\to\mathbb{R}$ be a cost function. The current section is about ergodic optimal control problems of the form

Here, $X$ denotes a martingale solution of the SDE (5) that depends implicitly on the optimization variable $\mu$.

Assumption 1 neither excludes that $\Omega$ is unbounded nor that $\Omega=\mathbb{R}^{n_{x}}$. On such unbounded domains, however, a weak Has’minskiǐ-Lyapunov function is needed.

This section discusses methods for analyzing the expected value of the stage cost $\ell$ at the optimal steady state of the FPK. As a first step towards this goal, the following lemma exploits a variant of Has’minskiǐ’s theorem in order to guarantee the existence of a potentially suboptimal steady state and to derive an associated upper bound on the expected value of $\ell$.

Proof. The construction of $\rho_{\mathrm{s}}$ and $\mu$ in Lemma 1 satisfies

As such, the closed and convex cone $\mathcal{K}$ is compatible with Slater’s constraint qualification, recalling that $\rho_{\mathrm{s}}>0$ follows from the minimum principle for elliptic operators [15]. Consequently, feasibility and boundedness of (32) follows immediately from Lemma 1, as $(\rho,\nu)=(\rho_{\mathrm{s}},\mu\rho_{\mathrm{s}})$ is a feasible point with bounded objective value. In summary, (32) has the following properties:

1. 1.

   Slater’s (weak) constraint qualification holds,

2. 2.

   the objective function $\mathcal{L}$ is a convex perspective function that is—due to Assumption 1—strongly convex in the control $\nu$, and

3. 3.

   the operator $\mathcal{A}$ is uniformly elliptic.

Let $(\rho_{\infty},\nu_{\infty})$ denote the optimal steady state, as defined in Corollary 1, and let $\mu_{\infty}=\nu_{\infty}/\rho_{\infty}$ denote an associated optimal control law. It can be used to study the controlled parabolic FPK

on a possibly infinite time horizon $\mathbb{T}$. Here, $\rho_{0}\in\mathcal{P}_{0}(\rho_{\infty}^{-1})$ denotes an initial probability density that is assumed to be an element of the set

In this setting, the key idea for analyzing (40) is to introduce the energy functional

Its weak time derivative is given by

Note that this equation is well-known in the functional analysis and FPK literature [3, 15], where $E$ is often interpreted as energy (or entropy) that is dissipated during the evolution of the FPK. Alternatively, by translating this to control engineering language, $E$ could also be called a Lyapunov function for (40).

Proof. FPKs have been analyzed by so many authors that is impossible to list them all. As such, [3] and [15] should be regarded as selected examples for two out of a long list of articles in which proofs of variants of the current lemma can be found. Nevertheless, because this is important for the developments in this paper, it appears appropriate to recall the following three main arguments why the above statement holds.

1. 1.

   First, if $\Omega_{\mathbb{T}}$ is bounded, a unique weak solution $\rho\in\mathcal{P}(\rho_{\infty}^{-1})$ of (40) exists, since it is a uniformly parabolic linear PDE. Moreover, it follows from Harnack’s inequality that

   $\displaystyle\phi\ \overset{\mathrm{def}}{=}\ \frac{\rho}{\rho_{\infty}}\ \in\ H^{1}(\rho_{\infty})\qquad\text{on}\ \ \mathbb{T}.$

   (42)

   The details of this argument can be found in [15, Chapter 7], where the square-integrability of logarithmic gradients of FPKs is discussed in all detail.

2. 2.

   Next, (41) can be verified by direct computation,

   	$\displaystyle\frac{\mathrm{d}E}{\mathrm{d}t}$	$\displaystyle\ =\$	$\displaystyle 2\left\langle\dot{\rho},\phi-1\right\rangle\ =\ 2\left\langle\dot{\rho},\phi-1\right\rangle-\left\langle\dot{\rho}_{\infty},\phi^{2}\right\rangle$		(43)
   		$\displaystyle\overset{\eqref{eq::FPinf}}{=}$	$\displaystyle 2\int_{\Omega}\left(\rho F[\mu_{\infty}]^{\intercal}\nabla\phi-\epsilon\nabla\rho^{\intercal}\nabla\phi\right)\,\mathrm{d}x$
   			$\displaystyle-2\int_{\Omega}\left[\rho_{\infty}F[\mu_{\infty}]^{\intercal}\nabla\phi-\epsilon\nabla\rho_{\infty}^{\intercal}\nabla\phi\right]\phi\,\mathrm{d}x$
   		$\displaystyle\overset{\eqref{eq::phi}}{=}$	$\displaystyle-2\epsilon\int_{\Omega}\left\|\nabla\phi\right\|_{2}^{2}\rho_{\infty}\,\mathrm{d}x\ \leq\ 0,$

   where the equation in the second line follows by using that the steady state $\rho_{\infty}$ also satisfies (40) but with $\rho_{0}=\rho_{\infty}$ and $\dot{\rho}_{\infty}=0$. Thus, Gronwall’s lemma yields the a-priori energy estimate

   $\displaystyle\|\rho(t)-\rho_{\infty}\|_{L^{2}(\rho_{\infty}^{-1})}\ \leq\ \|\rho_{0}-\rho_{\infty}\|_{L^{2}(\rho_{\infty}^{-1})}\quad\text{on}\ \ \mathbb{T}.$

   (44)

3. 3.

   For the case that $\Omega$ and $\mathbb{T}$ are unbounded, one can proceed as in the proof of Lemma 1. This means that the above two steps are repeated in order to construct a sequence of solutions $\rho_{k}\in\mathcal{P}(\rho_{\infty}^{-1})$ of (40) on the bounded domains $\Omega_{\alpha_{k}}\times(0,T_{k})$ with increasing sequences $\alpha_{k}\to\infty$ and $T_{k}\to\infty$. Because the energy estimate (44) yields a uniform upper bound, the sequence $(\rho_{k})_{k\in\mathbb{N}}$ has a weakly convergent subsequence in the Hilbert space $\mathcal{P}(\rho_{\infty}^{-1})$, whose limit must be a weak solution of (40) and, by continuity, (41) holds, too.

In order to analyze the ergodic optimal control problem (14), an additional assumption is introduced.

This assumption is equivalent to requiring that $\ell[\mu_{\infty}]$ has a variance at the steady state $\rho_{\infty}$; that is,

If Assumption 1 and 2 hold and if $\Omega$ is bounded, this is clearly the case, and Assumption 3 automatically holds. Otherwise, if $\Omega$ is unbounded, such a variance exists under very mild growth conditions on $\ell$, as discussed in all detail in [15, Chapter 3], where upper bound estimates for $\rho_{\infty}$ can be found. Next, let

denote the value function of the stochastic control problem on the finite time horizon $\mathbb{T}=(0,T)$ with parameters $T>0$ and $\rho_{0}\in\mathcal{P}_{0}(\rho_{\infty}^{-1})$. Here, the optimization variables $\rho$ and $\nu$ are elements of the convex cone

By construction, Problem (14) is equivalent to computing the limit

The following lemma proves that $J^{\star}=\ell_{\infty}$, which, in turn, implies that $J^{\star}$ does not depend on $\rho_{0}$.

Proof. This proof is divided into two parts, where the first part verifies the claim that (45) has a unique minimizer while the second part is about the ergodic limit.