Control of Fractional Diffusion Problems via Dynamic Programming Equations

Abstract

We explore the approximation of feedback control of integro-differential equations containing a fractional Laplacian term. To obtain feedback control for the state variable of this nonlocal equation, we use the Hamilton–Jacobi–Bellman equation. It is well known that this approach suffers from the curse of dimensionality, and to mitigate this problem we couple semi-Lagrangian schemes for the discretization of the dynamic programming principle with the use of Shepard approximation. This coupling enables approximation of high-dimensional problems. Numerical convergence toward the solution of the continuous problem is provided together with linear and nonlinear examples. The robustness of the method with respect to disturbances of the system is illustrated by comparisons with an open-loop control approach.

Appendix A: Derivation of example with analytic solution

Let $q$ and $\tilde{b}$ be the exact solution and corresponding right-hand side of the fractional Poisson problem

$$\begin{aligned} \left\{ \begin{array}{rcll} (-\Delta )^{s}q({\varvec{\xi }}) &{}=&{} \tilde{b}({\varvec{\xi }}) &{}{\varvec{\xi }}\in D, \\ q({\varvec{\xi }})&{}=&{}0 &{}{\varvec{\xi }}\in D^{c}, \\ \end{array} \right. \end{aligned}$$

with $\Vert q\Vert _{L^{2}(D)}=1$. An example of such a pair when $D=(-1,1)$ is given by

$$\begin{aligned} \begin{aligned}&\tilde{b}({\varvec{\xi }}) = 2^{2s}\Gamma \left( 1+s\right) \frac{\Gamma (s+1/2)}{\Gamma (1/2)} \sqrt{\frac{\Gamma (2s+3/2)}{\Gamma (2s+1)\Gamma (1/2)}}, \\ {}&q({\varvec{\xi }}) = \sqrt{\frac{\Gamma (2s+3/2)}{\Gamma (2s+1)\Gamma (1/2)}} \left( 1-{\varvec{\xi }}^{2}\right) ^{s}_{+}. \end{aligned} \end{aligned}$$

Now, let $\phi ,\kappa :(0,T)\rightarrow \mathbb {R}$ such that $\phi (0)=1$, $\kappa (T)=0$. Let $U:=[a,b]$ with $a<0<b$ and let $\gamma ,\lambda >0$. Consider the following functions that will be used in the construction of the cost functional and the state equation

$$\begin{aligned} \begin{aligned}&y_{d}({\varvec{\xi }},t) :=\phi (t)q({\varvec{\xi }}) - \gamma \kappa '(t)q({\varvec{\xi }}) + \gamma \kappa (t) \tilde{b}({\varvec{\xi }}) + \lambda \gamma \kappa (t)q({\varvec{\xi }}), \\ {}&u_{d}(t) := {\text {proj}}_{U}(\kappa (t)) ,\\ {}&b({\varvec{\xi }},t) := \phi '(t)q({\varvec{\xi }}) + \phi (t)\tilde{b}({\varvec{\xi }}) - u_{d}(t)q({\varvec{\xi }}). \end{aligned} \end{aligned}$$

We first aim to minimize the finite time horizon cost functional

$$\begin{aligned} \begin{aligned} \mathcal {J}_{q}^{T}(y, u)&:= \frac{1}{2} \int _0^{T} (\Vert y(\cdot ,\eta )-y_{d}(\cdot ,\eta )\Vert _{L^{2}(D)}^{2} + \gamma \Vert u(\cdot , \eta )\Vert _{L^{2}(D)}^2) e^{-\lambda \eta } d\eta \\ {}&= \frac{1}{2} \Vert y-y_{d}\Vert _{L_{\nu }^{2}(0,T;D)}^{2} + \frac{\gamma }{2} \Vert u\Vert _{L_{\nu }^{2}(0,T;D)}^{2}, \end{aligned} \end{aligned}$$

where we have set $\nu (t):=e^{-\lambda t}$, subject to

$$\begin{aligned} \left\{ \begin{array}{rll} \partial _t y({\varvec{\xi }},t) + (-\Delta )^{s}y({\varvec{\xi }},t) &{}= b({\varvec{\xi }},t) + u(t)q({\varvec{\xi }}) &{}({\varvec{\xi }},t) \in D\times (0,T){,} \\ y({\varvec{\xi }},t)&{}=0 &{}({\varvec{\xi }},t)\in D^{c} \times (0,T) , \\ y({\varvec{\xi }},0)&{}=q({\varvec{\xi }}) &{}{\varvec{\xi }}\in D. \end{array} \right. \end{aligned}$$

(42)

For fixed initial condition $q$, we can write the first order optimality conditions. Let $j(u):=\mathcal {J}_{q}^{T}(\mathcal {S}u,u)$, where $\mathcal {S}$ is the solution operator of the state equation Eq. (42) above. The operator $\mathcal {S}$ is also called the fractional control-to-state operator

$$\mathcal {S}: L^{2}(0,T;D) \rightarrow \mathbb {V}$$

and $\mathcal {S}u=z(u)$, where z(u) solves Eq. (42).

Then, the optimal control must satisfy the variational inequality

$$\begin{aligned} \overline{u}={\text {argmin}} j(u) \Leftrightarrow ( j'(\overline{u}), u-\overline{u})\ge 0 \quad \forall u \in \mathcal {U}_{ad}. \end{aligned}$$

The inequality can be written in the equivalent form

$$\begin{aligned} (\mathcal {S}^{*}(\mathcal {S}\kern0.1500em\overline{u}-y_{d}) + \gamma q\overline{u},qu-q\overline{u})_{L_{\nu }^{2}((0,T);D)}\ge 0, \end{aligned}$$

where $\mathcal {S}^{*}$ is the adjoint solution operator. Hence, $\overline{p}:=\mathcal {S}^{*}(\mathcal {S}\kern0.1500em\overline{u}-y_{d})$ is the solution to

$$\begin{aligned} \left\{ \begin{array}{rll} -\partial _t \overline{p}({\varvec{\xi }},t) +\lambda \overline{p} + (-\Delta )^{s}\overline{p}({\varvec{\xi }},t) &{}= \overline{y}({\varvec{\xi }},t)-y_{d}({\varvec{\xi }},t) &{}({\varvec{\xi }},t) \in D\times (0,T), \\ \overline{p}({\varvec{\xi }},t)&{}=0 &{}({\varvec{\xi }},t)\in D^{c} \times (0,T), \\ \overline{p}({\varvec{\xi }},T)&{}=0 &{}\text {for all }{\varvec{\xi }}\in D. \end{array} \right. \end{aligned}$$

with $\overline{z} = \mathcal{S}\kern0.1500em\overline{u}$. Then, the following inequality

$$\begin{aligned} 0\le (\overline{p} + \gamma q\overline{u},qu-q\overline{u})_{L_{\nu }^{2}(0,T;D)} = ((\overline{p},q)_{L^{2}(D)}q + \gamma q\overline{u},qu-q\overline{u})_{L_{\nu }^{2}(0,T;D)} \end{aligned}$$

implies that $\overline{u} = {\text {proj}}_{U}\left( -\frac{1}{\gamma }(\overline{p},q)_{L^{2}(D)}\right)$. If we set

$$\begin{aligned} y^{*}({\varvec{\xi }},t):=\phi (t)q({\varvec{\xi }}), \quad p^{*}({\varvec{\xi }},t):=-\gamma \kappa (t)q({\varvec{\xi }}), \quad u^{*}(t):=u_{d}(t), \end{aligned}$$

we obtain

$$\begin{aligned} \begin{aligned} y^{*}({\varvec{\xi }},0)&= q({\varvec{\xi }}), \\ \partial _t y^{*}({\varvec{\xi }},t) + (-\Delta )^{s}y^{*}({\varvec{\xi }},t)&= \phi '(t)q({\varvec{\xi }}) + \phi (t)\tilde{b}({\varvec{\xi }}) = b({\varvec{\xi }},t)+u^{*}({\varvec{\xi }},t), \end{aligned} \end{aligned}$$

where we have used that q is the solution of the fractional Poisson equation. Hence, $y^{*}$ is the state corresponding to the control $u^{*}$ and the initial condition $q$. Moreover,

$$\begin{aligned} \begin{aligned} p^{*}({\varvec{\xi }},T)&=0,\\ -\partial _t p^{*}({\varvec{\xi }},t) + \lambda p^{*}({\varvec{\xi }},t) + (-\Delta )^{s}p^{*}({\varvec{\xi }},t)&= \gamma \kappa '(t)q({\varvec{\xi }}) -\lambda \gamma \kappa (t)q({\varvec{\xi }}) -\gamma \kappa (t)\tilde{b}({\varvec{\xi }}) \\&= y^{*}({\varvec{\xi }},t)-y_{d}({\varvec{\xi }},t). \end{aligned} \end{aligned}$$

and hence $p^{*}$ solves the adjoint equation with right-hand side $y^{*}-y_{d}$. Finally, ${\text {proj}}_{U}\left( -\frac{1}{\gamma }(p^{*},q)_{L^{2}}\right) = {\text {proj}}_{U}\left( \kappa (t)\right) = u^{*}(t).$ Therefore, for the initial condition $q$, the optimal control is $u^{*}$, and the optimal state is $y^{*}$.

Let us now link the problem to the infinite time horizon framework. In order to do so, let us choose $T_{0}>0$ and $\kappa$ such that $\kappa (t)=0$ for $t\ge T_{0}$. For $T>T_{0}$, the previous construction gives the solution of the optimal control problem on $(0,T)$. On the other hand, we have

$$\begin{aligned} \begin{aligned} \mathcal {J}_{q}^{\infty }(y^{*},u^{*})&= \mathcal {J}_{q}^{T}(y^{*},u^{*}) + \frac{1}{2} \Vert y^{*}-y_{d}\Vert _{L_{\nu }^{2}(T,\infty ;D)}^{2} + \frac{\gamma }{2} \Vert u^{*}\Vert _{L_{\nu }^{2}(T,\infty ;D)}^{2} \\&= \mathcal {J}_{q}^{T}(y^{*},u^{*}) + \frac{\gamma ^{2}}{2} \left[ \Vert \kappa '\Vert _{L_{\nu }^{2}(T,\infty )}^{2} + \Vert \kappa \Vert _{L_{\nu }^{2}(T,\infty )}^{2}\Vert \tilde{b}\Vert _{L^{2}(D)}^{2}\right. -\lambda (\kappa ',\kappa )_{L_{\nu }^{2}(T,\infty )}\\&-\left. 2(\kappa ',\kappa )_{L_{\nu }^{2}(T,\infty )} (q,\tilde{b})_{L^{2}(D)} +2\frac{\lambda }{\gamma }\Vert \kappa \Vert _{L_{\nu }^{2}(T,\infty )}^{2}(q,\tilde{b})_{L^{2}(D)} + \lambda ^2 \Vert \kappa \Vert _{L_{\nu }^{2}(T,\infty )}\right] \\ {}&\qquad + \frac{\gamma }{2} \Vert {\text {proj}}_{U}(\kappa )\Vert _{L_{\nu }^{2}(T,\infty )}^{2} \\&= \mathcal {J}_{q}^{T}(y^{*},u^{*}), \end{aligned} \end{aligned}$$

where we have used that $\kappa$ and $\kappa '$ are zero on $(T,\infty )$. Therefore,

$$\begin{aligned} \min \mathcal {J}_{q}^{\infty }(y,u) \le \min \mathcal {J}_{q}^{T}(y,u), \end{aligned}$$

but the inverse inequality also holds, since $\mathcal {J}_{q}^{T}(y,u)\le \mathcal {J}_{q}^{\infty }(y,u)$. Hence, the pair $(y^{*},u^{*})$ is also optimal for the infinite time horizon case.