Data-driven initialization of deep learning solvers for Hamilton-Jacobi-Bellman PDEsfootnoteinfo

[footnoteinfo]DK was supported
by the UK Engineering and Physical Sciences Research Council (EPSRC) grants
EP/V04771X/1, EP/T024429/1, and EP/V025899/1.

First]A. Borovykh
Second]D. Kalise
Third]A. Laignelet
Third]P. Parpas

[First]Warwick Business School, University of Warwick, Coventry CV4 7AL, UK (e-mail<EMAIL_ADDRESS>[Second]Department of Mathematics, Imperial College London,
London SW7 2AZ, UK (e-mail<EMAIL_ADDRESS>[Third]Department of Computing, Imperial College London,
London SW7 2AZ, UK (e-mail: {alexis.laignelet18<EMAIL_ADDRESS>
A deep learning approach for the approximation of the Hamilton-Jacobi-Bellman partial differential equation (HJB PDE) associated to the Nonlinear Quadratic Regulator (NLQR) problem. A state-dependent Riccati equation control law is first used to generate a gradient-augmented synthetic dataset for supervised learning. The resulting model becomes a warm start for the minimization of a loss function based on the residual of the HJB PDE. The combination of supervised learning and residual minimization avoids spurious solutions and mitigate the data inefficiency of a supervised learning-only approach. Numerical tests validate the different advantages of the proposed methodology.

Hamilton-Jacobi-Bellman PDE, NLQR, supervised learning, residual minimization.

§ INTRODUCTION

Over the last years there has been remarkable progress on the solution of high-dimensional nonlinear optimal control problems, see e.g. [Dolgov et al., 2021, Onken et al., 2021, Kunisch, Karl and Walter, Daniel, 2021, Meng et al., 2022], partly fuelled by the effectiveness of deep learning methods for solving nonlinear PDEs including HJB equations (see e.g. [rai, 2019, Sirignano and Spiliopoulos, 2018, Han et al., 2018, Güler et al., 2019]). In its general form, unsupervised deep learning-based PDE solvers aim at minimizing a loss function in which the interior and boundary constraints from the PDE are incorporated (see e.g. eq. (<ref>)).
Certain challenges can still arise; the works of e.g. [van der Meer et al., 2020, Wang et al., 2021] have shown that the accuracy of the solution obtained by the neural network heavily relies on the weighting between the different terms (interior, boundary) in the loss function. As also noted in [van der Meer et al., 2020] a long training time can be required and the number of domain points one has to sample can become large. Furthermore, if the PDE solutions are non-unique, using solely the residual can lead to an incorrect solution; we show this for the linear HJB equation in Fig. <ref>. As an alternative, supervised learning methods can be employed (see e.g. [Azmi et al., 2021, Liu et al., 2019, Nakamura-Zimmerer et al., 2021]) where an optimal control solver is used to generate the solution dataset and a machine learning model is then trained on this set of data in a supervised manner. Alternatively, in [Albi et al., 2021] the solutions of the state-dependent Riccati equation (SDRE) (see e.g. [Banks et al., 2007]) are used as the initial dataset. A downside of such supervised settings is the propagation of errors from the initial solver into the model learned by the neural network (see Fig. <ref> where the discrepancy between the residual loss and the SDRE-generated solution is highlighted).

In this work we focus on the synthesis of feedback laws for optimal nonlinear stabilization problems, that is, by solving the HJB PDE associated to the NLQR problem. We propose a neural network-based approximation for the value function where the training is done using a two-step learning approach: in the first step (i) the parameters of the neural network are initialized on the state-value dataset obtained from solving the SDREs, and in the second step (ii) the residual coming from the PDE constraints is minimized. The benefits of this approach include:

* by first pre-training on the data, our method ensures the convergence to the correct solution; this is particularly relevant in a setting where the solution to the HJB PDE is non-unique (see e.g. Figure <ref>),

* by using the residual-minimization in the second learning step our method allows to mitigate the errors in the SDRE solution and allows to reach higher-accuracy solutions than data-only (see e.g. Figure <ref>),

* the two-step learning approach decreases the amount of data points required to reach a satisfactory accuracy compared to training on the dataset only (see e.g. Figure <ref>).

§ METHODOLOGY

§.§ Infinite horizon nonlinear optimal stabilization

Let $\by(t)\in\mathbb{R}^n$  and $\bu(t)\in\mathbb{R}^m$ be the state and control signal, respectively, of a system represented in control-affine form

\begin{align}\label{eq:dynamics_y}
\frac{d}{dt}\mathbf{y}(t)= \mathbf{f}(\by(t))+\mathbf{g}(\by(t))\bu(t), \;\; \by(0) = \bx\,.
\end{align}

We consider the following infinite horizon optimal control problem:

\begin{align}\label{eq:objective_u}
\min_{\bu(\cdot)\in\mathcal{U}}J(\bu;\bx):=\frac12\int_{0}^{\infty}\by^{\top}Q\by+\bu^{\top}R\bu \,dt\,,
\end{align}

subject to (<ref>), where the dynamics $\bff:\mathbb{R}^n\rightarrow\mathbb{R}^n$ and $\bg:\mathbb{R}^n\rightarrow\mathbb{R}^{n\times m}$ are assumed to be continuously differentiable, and the running cost matrices are given by $Q\in\mathbb{R}^{n\times n}, Q\succeq 0\,,$ and  $R\in\mathbb{R}^{m\times m}, R\succ 0$. We set $\mathcal{U}\equiv L^2([0,\infty);\R^m)$.

The formulation above corresponds to the NLQR problem, which arises in the synthesis of optimal feedback laws for the stabilization of nonlinear dynamics. The solution of this problem follows standard dynamic programming argument, see e.g., [Bardi et al., 1997]. Defining the value function of the problem

\begin{align}\label{eq:optimal_cost}
V(\mathbf{x})=\inf_{\mathbf{u}(\cdot)} J(\mathbf{u};\mathbf{x})\,,
\end{align}

the optimal control is expressed in feedback form as

\begin{align}\label{eq:opt_control}
\mathbf{u}^*(\mathbf{x})=-R^{-1}\mathbf{g}(\mathbf{x})^{\top}\nabla V(\mathbf{x})\,,
\end{align}

where $V(\bx)$ satisfies the Hamilton-Jacobi-Bellman PDE

\begin{align}\label{eq:hjb_pde}
& \mathcal{N}(\mathbf{x},V):=-\frac{1}{2} \nabla V(\mathbf{x})^{\top}\mathbf{g}(\mathbf{x})R^{-1}\mathbf{g}^{\top}(\mathbf{x})\nabla V(\mathbf{x})\\
& \qquad + \nabla V(\mathbf{x})^{\top}\mathbf{f}(\mathbf{x})+\frac12\mathbf{x}^{\top}Q\mathbf{x}=0\,,
\end{aligned}
\end{align}

in addition to the boundary condition $V(\mathbf{0})=0$.

§.§ Neural network model

In this work, we are interested in approximating the value function using deep neural networks. For this, we consider a fully connected feedforward neural network given by

\begin{align}
\hat V (\mathbf{x};\theta) = f_L\circ f_{L-1}\circ\cdots\circ f_1(\mathbf{x})
\end{align}

where $\mathbf{x}$ is the model input, $L$ is the number of layers and $f_i(\mathbf{x})=\sigma(W^i\mathbf{x}+\mathbf{b}^i)$, where $W^i$ and $b^i$ are matrices and vectors of prescribed dimensions, and $\sigma$ is a nonlinear activation function applied component-wise. The unknown parameters of the network $\{W^i,b^i\}_{i=1}^L$ are grouped in the vector $\theta$. The neural network is trained upon a set of collocation points $\{\bx^i\}_{i=1}^N$ sampled uniformly over a domain of interest. The training of the neural network in based on the two-step learning procedure presented in Section <ref>.

§.§ Synthetic data generation for supervised learning

In our two-step learning procedure $\hat V (\mathbf{x};\theta)$ is first trained following a supervised learning gradient-augmented approach, that is, upon a collection of values
\[\{\bx^i,V(\bx^i),\nabla V(\bx^i)\}_{i=1}^{N_1}\,,\]
approximating pointwise valuations of the true value function and its gradient. This dataset is generated synthetically by resorting to the link between the value function of the control problem and optimality conditions.  One possibility to generate the dataset in a finite horizon setting is to use Pontryagin's Maximum Principle to obtain first-order optimality conditions for (<ref>)-(<ref>) in the form of a two-point boundary value problem, separately for each initial condition $\bx^i$. However, in the infinite horizon optimal stabilization problem setting studied in this work, such an interpretation is not readily available, as the two-point boundary value problem has an asymptotic terminal condition for $t\to\infty$. We therefore rely on the use of a state-dependent Riccati equation approach [Banks et al., 2007] which provide a stabilizing feedback control through a sequential solution of the algebraic Riccati equations (ARE) along a trajectory.

§.§.§ State-Dependent Riccati Equation

For nonlinear dynamics, using the state-dependent Riccati equation (SDRE) approach is an effective method to generate a suboptimal, yet stabilizing solution which can be interpreted as an approximation of the HJB PDE (<ref>)   (see e.g. [Jones and Astolfi, 2020, Albi et al., 2021]). This approach is based on an extended linearization of the dynamics, expressed in semilinear form $f(\bx)=A(\bx)\bx$. Without loss of generality let $V(\bx)=\frac{1}{2}\mathbf{x}^{\top}P(\mathbf{x})\mathbf{x}$. Note that $\nabla V(\mathbf{x})=P(\bx)\mathbf{x} + \phi(\mathbf{x})$ where  $\phi\in\mathbb{R}^n$ is a `discrepancy' term arising from the differentiation of the non-linearity in $P(\mathbf{x})$ and its entries are given by

\begin{align*}
&\phi(\mathbf{x})_k = \sum_{i,j=1}^nx_ix_j\nabla_{x_k}P_{ij}(\mathbf{x})\,,\quad k=1,\ldots,n\,.
\end{align*}

In the SDRE approach, we approximate $\nabla V(\mathbf{x})\approx P(\mathbf{x})\bx$ so that the HJB PDE for the NLQR simplifies to

\begin{align}\label{eq:sdre}
\begin{aligned}
& -\frac{1}{2}P(\mathbf{x})g(\mathbf{x})R^{-1}g(\mathbf{x})^{\top}P(\mathbf{x})+P(\mathbf{x})A(\mathbf{x}) \\
& \qquad +A(\mathbf{x})^{\top}P(\mathbf{x})+Q=0,
\end{aligned}
\end{align}

For a fixed $\mathbf{x}$, the above equation amounts to solving an ARE. This suggests an implementation in a model predictive control fashion. Specifically, given a current state $\mathbf{\bar x}$ we solve (<ref>) for $P(\mathbf{\bar x})$ by fixing every operator. Then using (<ref>) gives
$\mathbf{u}(\mathbf{x})=-g^{\top}(\mathbf{\bar x})P(\mathbf{\bar x})\mathbf{x}$ and consequently using this in the dynamics of (<ref>) gives the next value of the state. Similarly as in the ARE setting, along a trajectory we generate data for $V(\bx)=\frac12 \bx^{\top}P(\bx)\bx$ and $\nabla V(\bx)=P(\bx)\bx$. However, as long as $\phi(\bx)\neq 0$, this is a suboptimal approximation to the solution of the HJB PDE. Through supervised learning, we train a neural network which is subsequently used as a warm start of the training of a neural network which will directly approximate the solution of the HJB PDE.

§.§ The training methodology

We propose a two-step learning approach. A first neural network is trained using supervised learning with a synthetic dataset generated from the SDRE approach. In a second step, a neural network, initialized with the parameters of the first step, is trained solely for minimizing the residual of the HJB PDE (<ref>).  The idea behind such a solver is that training first on the synthetic dataset allows the residual minimization to start with a set of weights closer to the basin of attraction of the correct solution. The residual minimisation is then necessary to reduce the approximation errors induced by the data generation from the SDRE approach.

§.§.§ Data-driven initialization

The SDRE approach generates a gradient-augmented dataset $\{\bx^i,V(\bx^i),\nabla V(\bx^i)\}_{i=1}^{N_1}$. Note that the gradient information is generated at virtually no extra computational cost. Consider the following loss function,

\begin{align}\label{eq:data_loss}
\begin{aligned}
\mathcal{L}^{\textnormal{dat}}(\theta) = &\lambda_1\frac{1}{N_1} \sum_{i=1}^{N_1}||V(\bx^i)-\hat V(\bx^i;\theta)||_2^2 \\
&+ \lambda_2\frac{1}{N_1} \sum_{i=1}^{N_1}||\nabla V(\bx^i)-\nabla \hat V(\bx^i;\theta)||_2^2,
\end{aligned}
\end{align}

where $N_1$ is the number of data points generated using the SDRE solver and $\lambda_1$ and $\lambda_2$ are the weights associated to the different loss function terms.
Our first step then amounts to initialising the neural network parameters by finding,

\begin{align}
\theta^{\textnormal{dat}} = \arg\min_{\theta} \mathcal{L}^{\textnormal{dat}}(\theta).
\end{align}

We remark that $\nabla \hat V(\bx;\theta)$ is obtained by differentiating the neural network with respect to the input $\mathbf{x}$.

§.§ Residual loss function

Based on (<ref>), the residual loss function is given by,

\begin{align}\label{eq:residual_loss}
\mathcal{L}^{\textnormal{res}}(\theta) = \frac{1}{N_2}\sum_{i=1}^{N_2}||\mathcal{N}(\mathbf{x}^i,\hat V(\bx^i;\theta))||_2^2,
\end{align}

where $N_2$ is the number of collocation points $\mathbf{x}^1,...,\mathbf{x}^{N_2}$ which we sample on the interior of the domain.
As before, the gradient terms in this loss function are obtained by differentiating the neural network output $\hat V$ with respect to state variables.
After obtaining $\theta^{dat}$ from minimizing (<ref>), these parameters are used as the initialization of the training of a second neural network based on the residual minimization:

\begin{align}
\theta^{\textnormal{res\_valid}} = \arg\min_{\theta} \mathcal{L}^{\textnormal{res}}(\theta)\,.
\end{align}

§.§ Benefits of the two-step learning approach

We illustrate, using a toy example, the benefits of warm start for residual minimization. Consider the linear quadratic control problem with $A=I$, $B=I$, $Q=I$ and $R=0.2I$ with $n=2$. The solution $P$ of the corresponding ARE is a diagonal matrix with elements $\alpha_1 \approx \pm 0.69$ and/or $\alpha_2 \approx\pm -0.29$, see Figure <ref> for the different solutions $V(\bx)$. By the positive-definiteness assumption the only valid solution is $P=\alpha_1 I$. However, if a neural network was trained on the residual only, four global minima exist (the solutions for $P$ with $P$'s elements different combinations of the roots) and convergence towards the value function will depend on a correct initialization[Note that in dimension $n$, the number of possible solutions in this example grows to $2^n$.].

Using Xavier initialization [Glorot and Bengio, 2010] the initial weights are centred around zero and hence gradient descent is most likely to converge to the minimum at $(-0.29, -0.29)$ (see the left-hand side of Fig <ref>).
To ensure we converge to the valid solution, we can enrich the loss function with a data-driven term, similar to $\mathcal{L}^{\textnormal{dat}}(\theta)$; this term guides convergence to the correct solution at $(0.69, 0.69)$ (see the right-hand side of Fig. <ref>).

group style=
group name=my plots,
group size=2 by 4,
ylabels at=edge left,
horizontal sep=1.8cm

mesh/ordering=x varies,
] table[x index=0,y index=1, z index=2,col sep=comma] data/lqr_true_solutions.dat;

mesh/ordering=x varies,
] table[x index=0,y index=1, z index=3,col sep=comma] data/lqr_true_solutions.dat;

mesh/ordering=x varies,
] table[x index=0,y index=1, z index=5,col sep=comma] data/lqr_true_solutions.dat;

mesh/ordering=x varies,
] table[x index=0,y index=1, z index=4,col sep=comma] data/lqr_true_solutions.dat;

Possible solutions of the Hamilton-Jacobi equation in the LQR case with $A = B = Q = I$, and $R = 0.2I$.

This toy example motivates our combined approach of data and residual loss minimization. Using a residual-only approach inherits the multiplicity of (wrong) solutions of the HJB PDE associated to the NLQR problem. The downside of using both data and residual terms in a single loss function is the problem of choosing the right weights between the residual and data term ([van der Meer et al., 2020]).
Even more crucially, in non-linear settings the  approximate value function obtained from the SDRE dataset can significantly differ from the value function that satisfies the residual (see also Figure <ref>); in such a setting, minimizing the data term and the residual term at the same time is not possible. To ensure our framework is robust and efficient in nonlinear settings, we thus propose the two-step learning approach.

title=Residual only,
width=, height=,

mesh/ordering=x varies,

[opacity=1] graphics[xmin=-0.6,xmax=1,ymin=-0.6,ymax=1] data/residual_only.pdf;
[only marks, mark=*, mark size = 1, black]
[only marks, mark=*, mark size = 1.5, red]
[only marks, mark=x, mark size = 3, black]
[ draw,
mark=between positions 0.4 and 0.9 step 3mm with latex,
postaction=decorate, black] (0,0) to (-0.29,-0.29);
[ draw, black, dashed] (0.2,-0.6) to (0.2,1);
[ draw, black, dashed] (-0.6, 0.2) to (1,0.2);

title=Residual + data,
width=, height=,

mesh/ordering=x varies,

[opacity=0.8] graphics[xmin=-0.6,xmax=1,ymin=-0.6,ymax=1] data/residual_data.pdf;
[only marks, mark=*, mark size = 1.5, red]
[only marks, mark=x, mark size = 3, black]

[ draw,
mark=between positions 0.2 and 0.95 step 4mm with latex,
postaction=decorate, black] (0,0) to (0.69,0.69);

Loss functions in the case of residual-only and a linear combination of the residual and the data term. The global minimum for the valid value function is represented by a red dot. A loss function combining the data and residual has a unique solution, while a residual-only approach inherits the multiple, but wrong, solutions of the original HJB PDE.

§ EXPERIMENTS

§.§ Implementation details

§.§.§ Neural network configuration and loss function

The architecture we consider is a fully connected neural network with the sigmoid activation function $\sigma(x)=1/(1+e^{-x})$ with 3 hidden layers. For the two-dimensional problem we set 20 neurons per layer while for the 10-dimensional case we use 50 neurons per layer. The neural network is trained using the Adam optimizer [Kingma and Ba, 2014].

For the 2D examples, 20 SDRE data points are sampled for supervised learning and 50 collocation points are used to minimize the HJB PDE residual. The training procedure on the dataset is the following: we use $10^3$ iterations with a learning rate of $0.01$ followed by $10^{3}$ iterations with a learning rate of $10^{-3}$. When training simultaneously on the data and gradient we set $\lambda_1,\lambda_2=1$.
For the residual minimization, the training is first done for $2\times 10^3$ iterations with a learning rate of $0.01$, followed by $4\times 10^{3}$ iterations with a learning rate of $10^{-3}$.

For the Cucker-Smale model, 500  data points are sampled, and $5\times 10^3$ points are used for training on the residual. The pre-training uses $2\times 10^3$ iterations with  a learning rate of $0.01$ followed by $4\times 10^3$ iterations with $10^{-3}$. The residual minimization is done using $10^4$ iterations with  a learning rate of $10^{-3}$ followed by $2\times 10^4$ iterations with $10^{4}$.

§.§.§ Data generation

For the 2D examples, the ARE/SDRE is solved analytically to obtain a closed-form solution for $P(\mathbf{x})$ which is then used in the expressions for $V(\mathbf{x})$ and $\nabla V(\mathbf{x})$.

In the case of the Cucker Smale model, the SDRE is solved numerically using the SciPy Linear Algebra package and $\nabla V(\mathbf{x})=P(\mathbf{x})\mathbf{x}$ is used for the derivative. The collocation points are uniformly sampled in the space.

§.§.§ Evaluation

The evaluation of the accuracy of the final solution is done by computing the mean squared error (MSE) between the neural network on newly generated interior domain points and the exact solution. For the 2D tests, the exact solution is obtained either through a closed-form solution, whenever available, or using the numerical solution of a semi-Lagrangian scheme over a fine grid [Alla et al., 2015]. The final MSE's are computed as the median over 10 runs of the optimization.

§.§ A 2D linear example

We revisit the same linear quadratic toy example discussed in Section <ref>. The left-hand side of Fig. <ref> shows the MSE for three different training methods: i) training on value function data only, ii) training on data on the value function and the gradient, and iii) training on residual only. By minimizing over the residual only, the neural network converges to the incorrect solution (as was also shown in the toy example in Fig. <ref>); by incorporating the gradient of the value function into the loss, a more accurate solution can be obtained. In the left-hand side of Fig. <ref> it is shown that two-step learning requires just 10 data points to obtain a solution that is of accuracy $10^{-6}$, while when learning only over the data, one would require 100 data points to achieve the same error. Comparing the left-hand sides of Fig <ref> and  <ref> we also note the two-step learning can converge to a solution with a smaller MSE than the data-only training. Finally, the left-hand side of Table <ref> shows the MSE for the data-only (using value function and gradient) and two-step learning. It is clear that when using more than 5 data points a more accurate solution is obtained using two-step learning.

group style=
group name=my plots,
group size=2 by 1,
ylabels at=edge left

ylabel shift = -5 pt,

[title=2D linear,ylabel=MSE]
[color=graph_green, thick]
table[x index=0,y index=1,col sep=comma] data/linear_all_methods.dat;

[color=graph_red, thick]
table[x index=0,y index=2,col sep=comma] data/linear_all_methods.dat;
[color=graph_petrol, thick]
table[x index=0,y index=4,col sep=comma] data/linear_all_methods.dat;

[title=2D nonlinear,legend style=at=(0.8,-0.3), anchor=north east,legend columns=3]
[color=graph_green, thick]
table[x index=0,y index=1,col sep=comma] data/non_linear_all_methods.dat;

[color=graph_red, thick]
table[x index=0,y index=2,col sep=comma] data/non_linear_all_methods.dat;
[color=graph_petrol, thick]
table[x index=0,y index=4,col sep=comma] data/non_linear_all_methods.dat;

Mean squared error of different training methods for 2D linear and 2D non linear case (computed over $100\times 100$ points in the domain).

group style=
group name=my plots,
group size=2 by 1,
ylabels at=edge left

ylabel shift = -5 pt,

xticklabels=0, 5, 10, 20, 50, 100,
xlabel=Count data points

[title=2D linear,ylabel=MSE]
[mark=x, mark size=3pt, color=graph_green, thick]
table[x index=0,y index=2,col sep=comma] data/linear_results.dat;
[mark=x, mark size=3pt, color=graph_petrol, thick]
table[x index=0,y index=3,col sep=comma] data/linear_results.dat;

[title=2D nonlinear,legend style=at=(0.4,-0.3), anchor=north east,legend columns=2]
[mark=x, mark size=3pt, color=graph_green, thick]
table[x index=0,y index=2,col sep=comma] data/non_linear_results.dat;
[mark=x, mark size=3pt, color=graph_petrol, thick]
table[x index=0,y index=3,col sep=comma] data/non_linear_results.dat;

Residual (pre-trained)

Mean squared error between training on dataset and two step learning (computed over $100\times 100$ points in the domain).

Mean squared error for 2D cases comparing the supervised learning on the dataset and the two steps learning.

2c2D linear     2c2D non linear

(lr)2-3 (lr)4-5
Point         Data/grad       2 steps       Data/grad       2 steps
0     0.20288865     0.04377991     2.181217     1.176169

5     0.00598487    0.02723277      0.251505    1.179283

10    0.00043885     0.00000109    0.04561     0.754814

20     0.00004997     0.00000121     0.013527     0.000343

50     0.00000280     0.00000153     0.00826     0.000206

100     0.00000109     0.00000132     0.008112     0.000156

§.§ A 2D nonlinear example

We next consider the dynamics given by

\[\dot \bx=\begin{bmatrix}
0 & 1  \\
\epsilon x_1^2 & 0
\end{bmatrix}\bx+\begin{bmatrix}
0  \\
\end{bmatrix}u(t)\,,\]

with $R = Q = I$. Note that the coefficient $\epsilon$ controls the intensity of the nonlinear term. The solution of the SDRE in this case is given by:

\begin{align*}
\sqrt{\epsilon^2 x_1^4 + 1}\sqrt{1 + 2P_{12}} & P_{12} \\
P_{12}  & \sqrt{1 + 2P_{12}}
\end{bmatrix}
\end{align*}

with $P_{12} = \epsilon x_1^2 + \sqrt{\epsilon^2 x_1^4 + 1}$.
Remember that the exact SDRE solution is an approximation of the original HJB PDE. In Fig. <ref> we show the residual generated by the exact SDRE solution when inserted into the HJB PDE as a function of $\epsilon$ - this shows that the more non-linear the problem is, the more relevant the consequent minimization of the residual will be. This also shows why it is not possible to combine both (<ref>) and (<ref>) into a single objective: the SDRE loss term and the residual loss term lead to different solutions and minimizing both at the same time leads to a decrease in accuracy compared to the two-step learning approach.
The right-hand side of Fig. <ref> shows that using no data points and training on the residual only leads to inaccurate solutions, while training on the dataset only does not result in an error lower than $10^{-2}$ due to the SDRE solution being only an approximate solution. Combining data-driven initialisation with the residual minimization results in the most accurate solution as shown in the right-hand side of Fig. <ref>: training on 20 points first, and then minimizing the residual leads to a solution with a $10^{-4}$ accuracy. This level of accuracy is not obtained by increasing the amount of data points used in the supervised learning step. This highlights the necessity of the residual minimization to improve the accuracy of the final solution. The right-hand side of Table <ref> validates these results. Finally, in Fig. <ref> we present the results from training on the residual only (the left-hand side) and the two-step learning (the right-hand side). The residual-only training results in an invalid solution. Based on the value of the residual loss in (<ref>) (the bottom plot) one could conclude that correct convergence has been obtained (the residual loss is sufficiently small); however as shown by the MSE computed by comparing the exact solution with the neural network output the network has converged to the incorrect solution. The value of two-step learning is clear from the right-hand side: the residual and MSE are both sufficiently low so that the solution obtained is close to the true solution.

[xlabel=Non linearity parameter $\epsilon$, ylabel=Hamilton Jacobi residual, width=0.9, height= 0.6,

try min ticks=4,

ylabel near ticks,
xticklabels=$10^{-3}$, $10^{-2}$, $10^{-1}$, $1$, $10$, $100$]

[mark=x, mark size=4pt, color=black, thick]
table[x index=0,y index=2,col sep=comma] data/non_linear_discrepancy.dat;

Influence of the intensity of the non linearity induced by $\epsilon$.

[group style=
group name=my plots,
group size=2 by 1,
ylabels at=edge left,

title= Exact solution of HJB PDE,
mesh/ordering=x varies,
] table[x index=0,y index=1, z index=6,col sep=comma] data/function_all.dat;

group style=
group name=my plots,
group size=2 by 4,
ylabels at=edge left,
horizontal sep=1.2cm

title=Residual only,
mesh/ordering=x varies,
] table[x index=0,y index=1, z index=2,col sep=comma] data/function_all.dat;

title= Two-step learning,
mesh/ordering=x varies,
] table[x index=0,y index=1, z index=3,col sep=comma] data/function_all.dat;

3[only marks, mark size = 1.2] table[x index=0,y index=1, z index=2,col sep=comma] data/dataset_coordinates.dat;

[color=graph_green, thick]
table[x index=0,y index=1,col sep=comma] data/non_linear_residual_solutions.dat;
[color=graph_petrol, thick]
table[x index=0,y index=2,col sep=comma] data/non_linear_residual_solutions.dat;

[color=graph_green, thick]
table[x index=0,y index=3,col sep=comma] data/non_linear_residual_solutions.dat;
[color=graph_petrol, thick]
table[x index=0,y index=4,col sep=comma] data/non_linear_residual_solutions.dat;

Training the residual from a random initialisation (left) and starting from a pre-training on a dataset (right). From top to bottom: the initialisation, the solution after training, the true solution and the value of the residual loss in (<ref>) and the MSE.

§.§ Cucker-Smale model

We conclude with a high-dimensional nonlinear test consisting of agent-based consensus dynamics, known as the Cucker-Smale model. We consider 5 agents having states $\mathbf{x}_i=(y_i,v_i)\in\mathbb{R}^2$ in $[-3,3]^{10}\in\mathbb{R}^{5}\times\mathbb{R}^{5}$, with dynamics given in semilinear form by

\begin{align}
\mathbf{\dot{x}}=\begin{bmatrix}
O_5 & I_5\\
O_5& A(\mathbf{y})
\end{bmatrix} \mathbf{x} + \begin{bmatrix}
\end{bmatrix}\bu\,,
\end{align}

where $O_5$ is an $5\times 5$ matrix of zeros and $I_5$ is the $5\times 5$ identity matrix and where

\begin{align}
-\frac{1}{5}\sum_{k=1}^5 \frac{1}{1+||y_i-y_k||^2} & \text{ for } i=j\,,\\
\frac{1}{5}\frac{1}{1+||y_i-y_j||^2} & \text{ otherwise}\,.
\end{cases}
\end{align}

We set $Q = \frac{1}{5}I_{10}$ and $R=I_5$. Note that the resulting HJB PDE is 10-dimensional. In Fig. <ref> we present the results for the position, velocity and control for i) residual-only training, ii) data and gradient-only training and iii) two-step learning. As before, it can be seen that training on just the residual  does not lead to the correct control. However, by first training on a dataset of 500 points, and then minimizing the residual, we are able to obtain a smooth stabilizing control. We remark that for this particular problem, 500 points is not enough for the supervised learning task, as in this data-only training the obtained control is not sufficiently smooth and many more points would be required to get a satisfying control that generalises well on unseen data. This highlights the advantage of our combined approach.

group style=
group name=my plots,
group size=3 by 2,
ylabels at=edge left,
horizontal sep=0.7cm

enlarge x limits=false,
max space between ticks=1000pt,

[title=Residual only,ylabel=Velocity, ylabel shift = -5 pt, ymin=-7, ymax=3]
[evaluate=as using (-6)*100/(10-6)] in 6,...,10

table[x index=0,y index= ] data/pde_residual_invalid_trajectories.dat;

[title=Data/gradient,ylabel shift = -5 pt, ymin=-1.7, ymax=3]
[evaluate=as using (-6)*100/(10-6)] in 6,...,10

table[x index=0,y index= ] data/pde_data_trajectories.dat;

[title=2 steps learning,ymin=-1.7, ymax=3]
[evaluate=as using (-6)*100/(10-6)] in 6,...,10

table[x index=0,y index= ] data/pde_residual_trajectories.dat;

[ylabel=Control, ylabel shift = -5 pt, xlabel=Time, ymin=-7, ymax=1]
[evaluate=as using (-1)*100/(5-1)] in 1,...,5

table[x index=0,y index= ] data/pde_residual_invalid_control.dat;

[ylabel shift = -5 pt, xlabel=Time, ymin=-7, ymax=1]
[evaluate=as using (-1)*100/(5-1)] in 1,...,5

table[x index=0,y index= ] data/pde_data_control.dat;

[xlabel=Time, ymin=-7, ymax=1]
[evaluate=as using (-1)*100/(5-1)] in 1,...,5

table[x index=0,y index= ] data/pde_residual_control.dat;

Trajectories and velocities of a group of particles under the control found by the neural network after (left) dataset-only training on 500 points and (right) pre-training on 500 points & minimising the residual.

§.§.§ Conclusion and discussion

We have shown how a two-step learning approach can result in efficient and accurate solutions to the HJB PDE. In our setting, the output of the neural network was set to directly be the value function $V(\mathbf{x})$. However, as in the SDRE derivation, it is possible to assume that $V(\bx)=\frac{1}{2}\mathbf{x}^{\top}P(\mathbf{x})\mathbf{x}$ and adapt the neural network architecture accordingly. Depending on the problem this can result in an efficient solver; we will address the choice of architecture for HJB PDE solvers in future work. Our results show that in a setting where the residual is sufficiently low for multiple solutions, the neural network converges to the one closest to initialization, which can be analyzed in the context of implicit bias [Chizat and Bach, 2020].

[rai, 2019]
Physics-informed neural networks: A deep learning framework for
solving forward and inverse problems involving nonlinear partial differential
J. Comput. Phys., 378, 686–707.

[Albi et al., 2021]
Albi, G., Bicego, S., and Kalise, D. (2021).
Gradient-augmented supervised learning of optimal feedback laws using
state-dependent riccati equations.
IEEE Contr. Syst. Lett, 6, 836–841.

[Alla et al., 2015]
Alla, A., Falcone, M., and Kalise, D. (2015).
An efficient policy iteration algorithm for dynamic programming
SIAM J. Sci. Comput., 37(1), A181–A200.

[Azmi et al., 2021]
Azmi, B., Kalise, D., and Kunisch, K. (2021).
Optimal feedback law recovery by gradient-augmented sparse polynomial
J. Mach. Learn. Res., 22(48), 1–32.

[Banks et al., 2007]
Banks, H., Lewis, B., and Tran, H.T. (2007).
Nonlinear feedback controllers and compensators: a state-dependent
riccati equation approach.
Comput. Optim. Appl., 37(2), 177–218.

[Bardi et al., 1997]
Bardi, M., Dolcetta, I.C., et al. (1997).
Optimal control and viscosity solutions of
Hamilton-Jacobi-Bellman equations, volume 12.

[Chizat and Bach, 2020]
Chizat, L. and Bach, F. (2020).
Implicit bias of gradient descent for wide two-layer neural networks
trained with the logistic loss.

[Dolgov et al., 2021]
Dolgov, S., Kalise, D., and Kunisch, K.K. (2021).
Tensor decomposition methods for high-dimensional
hamilton–jacobi–bellman equations.
SIAM J. Sci. Comput., 43(3), A1625–A1650.

[Glorot and Bengio, 2010]
Glorot, X. and Bengio, Y. (2010).
Understanding the difficulty of training deep feedforward neural
In Proceedings of the Thirteenth International Conference on
Artificial Intelligence and Statistics, 249–256.

[Güler et al., 2019]
Güler, B., Laignelet, A., and Parpas, P. (2019).
Towards robust and stable deep learning algorithms for forward
backward stochastic differential equations.

[Han et al., 2018]
Han, J., Jentzen, A., and Weinan, E. (2018).
Solving high-dimensional partial differential equations using deep
Proc. Nat. Acad. Sci., 115(34), 8505–8510.

[Jones and Astolfi, 2020]
Jones, A. and Astolfi, A. (2020).
On the solution of optimal control problems using parameterized
state-dependent riccati equations.
In 2020 59th IEEE Conference on Decision and Control (CDC),

[Kingma and Ba, 2014]
Kingma, D.P. and Ba, J. (2014).
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980.

[Kunisch, Karl and Walter, Daniel, 2021]
Kunisch, Karl and Walter, Daniel (2021).
Semiglobal optimal feedback stabilization of autonomous systems via
deep neural network approximation*.
ESAIM: COCV, 27, 16.

[Liu et al., 2019]
Liu, S., Oosterlee, C.W., and Bohte, S.M. (2019).
Pricing options and computing implied volatilities using neural
Risks, 7(1), 16.

[Meng et al., 2022]
Meng, T., Zhang, Z., Darbon, J., and Karniadakis, G.E. (2022).
Sympocnet: Solving optimal control problems with applications to
high-dimensional multi-agent path planning problems.

[Nakamura-Zimmerer et al., 2021]
Nakamura-Zimmerer, T., Gong, Q., and Kang, W. (2021).
Adaptive deep learning for high-dimensional hamilton–jacobi–bellman
SIAM J. Sci. Comput., 43(2), A1221–A1247.

[Onken et al., 2021]
Onken, D., Nurbekyan, L., Li, X., Fung, S.W., Osher, S., and Ruthotto, L.
A neural network approach for real-time high-dimensional optimal

[Sirignano and Spiliopoulos, 2018]
Sirignano, J. and Spiliopoulos, K. (2018).
Dgm: A deep learning algorithm for solving partial differential
J. Comput. Phys., 375, 1339–1364.

[van der Meer et al., 2020]
van der Meer, R., Oosterlee, C., and Borovykh, A. (2020).
Optimally weighted loss functions for solving pdes with neural

[Wang et al., 2021]
Wang, S., Teng, Y., and Perdikaris, P. (2021).
Understanding and mitigating gradient flow pathologies in
physics-informed neural networks.
SIAM Journal on Scientific Computing, 43(5), A3055–A3081.


# Initial-boundary value problem

Oqila Muhiddinova  Institute of Mathematics, Uzbekistan Academy of Science
Institute of Mathematics, Uzbekistan Academy of Science, Tashkent, 81 Mirzo
Ulugbek str. 100170<EMAIL_ADDRESS>
# Initial-boundary value problem for a subdiffusion equation with the Caputo
derivative

Oqila Muhiddinova  Institute of Mathematics, Uzbekistan Academy of Science
Institute of Mathematics, Uzbekistan Academy of Science, Tashkent, 81 Mirzo
Ulugbek str. 100170<EMAIL_ADDRESS>
###### Abstract.

We investigate an initial-boundary value problem for a time-fractional
subdiffusion equation with the Caputo derivatives on $N$-dimensional torus by
the classical Fourier method. Since our solution is established on the
eigenfunction expansion of elliptic operator, the method proposed in this
article can be used to an arbitrary domain and an elliptic operator with
variable coefficients. It should be noted that the conditions for the
existence of a solution to the initial-boundary value problem found in the
article cannot be weakened, and the article provides a corresponding example.

AMS 2000 Mathematics Subject Classifications : Primary 35R11; Secondary 74S25.
Key words: Time-fractional subdiffusion equation, the Caputo derivatives,
initial-boundary value problem, Fourier method, Liouville spaces.

+-

##  1\. Main result

The fractional derivative in the sense of Caputo of order $0<\rho<1$ of the
function $h(t)$ defined on $[0,\infty)$ has the form (see, for example, [1],
p. 14)

$D_{t}^{\rho}h(t)=\frac{1}{\Gamma(1-\rho)}\int_{0}^{t}\frac{\frac{d}{d\tau}h(\tau)d\tau}{(t-\tau)^{\rho}},\quad
t>0,$

provided the right-hand side exists. Here $\Gamma(\sigma)$ is Euler’s gamma
function. If in this definition we interchange differentiation and fractional
integration, then we get the definition of the Riemann-Liouville derivative:

$\partial_{t}^{\rho}h(t)=\frac{1}{\Gamma(1-\beta)}\frac{d}{dt}\int_{0}^{t}\frac{h(\tau)d\tau}{(t-\tau)^{\rho}},\quad
t>0.$

Note that if $\rho=1$, then fractional derivatives coincides with the ordinary
classical derivative of the first order:
$\partial_{t}h(t)=D_{t}h(t)=\frac{d}{dt}h(t)$.

Let $\mathbb{T}^{N}$ be $N$ -dimensional torus:
$\mathbb{T}^{N}=(-\pi,\pi]^{N}$, $N\geq 1$. We define by $C(\mathbb{T}^{N})$
and $C^{2}(\mathbb{T}^{N})$ a class of $2\pi$-periodic on each variable
$x_{j}$ functions $v(x)$ from $C(\mathbb{T}^{N})$ and $C^{2}(\mathbb{T}^{N})$
correspondingly. Let $A$ stand for a positive operator, defined on
$C^{2}(\mathbb{T}^{N})$ and acting as $Av(x)=-\Delta v(x)$, where $\Delta$ is
the Laplace operator.

Let $\rho\in(0,1)$ be a constant number. Consider the initial-boundary value
problem

(1.1) $D_{t}^{\rho}u(x,t)+Au(x,t)=f(x,t),\quad x\in\mathbb{T}^{N},\quad
0<t\leq T,$

(1.2) $u(x,0)=\varphi(x),\quad x\in\mathbb{T}^{N},$

where $f$ and $\varphi$ are given continuous functions.

######  Definition 1.1.

A function $u(x,t)\in C(\mathbb{T}^{N}\times[0,T])$ with the properties
$D_{t}^{\rho}u(x,t),A(x,D)u(x,t)\in C(\mathbb{T}^{N}\times(0,T])$, and
satisfying the conditions of problem (1.1) - (1.2) is called the classical
solution of the initial-boundary value problem.

Initial-boundary value problem (1.1) - (1.2) for various elliptic operators
$A$ have been considered by a number of authors using different methods (see,
for example, handbook [2]). In the book of A.A. Kilbas et al. [3] (Chapter 6)
there is a survey of works published before 2006. The case of one spatial
variable $x\in\mathbb{R}$ and subdiffusion equation with the elliptical part
$u_{xx}$ were considered for example in the monograph of A. V. Pskhu [1]
(Chapter 4, see references thesein).

The main method used in this work is the Fourier method. As far as we know, in
all previous papers this method was used to prove the existence of a
generalized solution of initial-boundary value problems for subdiffusion
equations. So for example, in the paper of Yu. Luchko [4] the author
constructed solutions by the eigenfunction expansion in the case of $f=0$ and
discussed the unique existence of the generalized solution to problem (1.1) \-
(1.2) with the Caputo derivative. In an arbitrary $N$-dimensional domain
$\Omega$ initial-boundary value problems for subdiffusion equations (the
fractional part of the equation is a multi-term and initial conditions are
non-local) with the Caputo derivatives has been investigated by M. Ruzhansky
et al. [5]. The authors proved the existence and uniqueness of the generalized
solution to the problem by the Fourier method. The authors of the paper [6]
investigated both problem (1.1) - (1.2) and the corresponding backward problem
with an arbitrary elliptic operator of the second order. To prove the
existence and uniqueness of the generalized solution the authors applied the
Fourier method.

The authors of papers [7] and [8] used the Fourier method to construct a
classical solution of the subdiffusion equations with the Riemann-Liouville
derivative and various elliptic operators.

Let us introduce some concepts and formulate the main result of the work.

Let $\hat{A}$ stand for the closure of operator $A$ in
$L_{2}(\mathbb{T}^{N})$. Then $\hat{A}$ is selfadjoint and it has a complete
(in $L_{2}(\mathbb{T}^{N})$) set of eigenfunctions $\\{\gamma e^{inx}\\}$,
$\gamma=\gamma(N)=(2\pi)^{-N/2}$, $n\in\mathbb{Z}^{N}$ and corresponding
eigenvalues $|n|^{2}=n_{1}^{2}+n_{2}^{2}+...+n_{N}^{2}$. Therefore, by virtue
of J. von Niemann theorem, for any $\tau>0$ one can introduce the power of
operator $\hat{A}$ as
$\hat{A}^{\tau}g(x)=\sum\limits_{n\in\mathbb{Z}^{N}}|n|^{\tau}g_{n}e^{inx}$,
where $g_{n}$ are Fourier coefficients:

$g_{n}=(2\pi)^{-N}\int\limits_{\mathbb{T}^{N}}g(x)e^{-inx}dx.$

The domain of definition of this operator is defined from the condition
$\hat{A}^{\tau}g(x)\in L_{2}(\mathbb{T}^{N})$ and has the form

$D(\hat{A}^{\tau})=\\{g\in
L_{2}(\mathbb{T}^{N}):\sum\limits_{n\in\mathbb{Z}^{N}}|n|^{2\tau}|g_{n}|^{2}<\infty\\}.$

If we denote by $L_{2}^{a}(\mathbb{T}^{N})$, $a>0$, the Liouville space with
the norm

(1.3)
$||g||^{2}_{L_{2}^{a}(\mathbb{T}^{N})}=\big{|}\big{|}\sum\limits_{n\in\mathbb{Z}^{N}}(1+|n|^{2})^{\frac{a}{2}}g_{n}e^{inx}\big{|}\big{|}^{2}_{L_{2}(\mathbb{T}^{N})}=\sum\limits_{n\in\mathbb{Z}^{N}}(1+|n|^{2})^{a}|g_{n}|^{2},$

then it is not hard to see that $D(\hat{A}^{\tau})=L_{2}^{\tau
m}(\mathbb{T}^{N})$.

Let $E_{\rho,\mu}$ be the two-parametric Mittag-Leffler function:

$E_{\rho,\mu}(t)=\sum\limits_{k=0}^{\infty}\frac{t^{k}}{\Gamma(\rho k+\mu)}.$

Here is the main result.

######  Theorem 1.2.

Let $a>\frac{N}{2}$ and $\varphi\in L^{a}_{2}(\mathbb{T}^{N})$. Moreover, let
$f(x,t)\in L^{a}_{2}(\mathbb{T}^{N})$ for $t\in[0,T]$ and
$||f(\cdot,t)||^{2}_{L_{2}^{a}(\mathbb{T}^{N})}\in C[0,T]$. Then there exists
a solution of initial-boundary value problem (1.1) - (1.2) and it has the form

(1.4)
$u(x,t)=\sum\limits_{n\in\mathbb{Z}^{N}}\bigg{[}\varphi_{n}E_{\rho,1}(-|n|^{2}t^{\rho})+\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})d\xi\bigg{]}e^{inx},$

which absolutely and uniformly converges on $x\in\mathbb{T}^{N}$ and for each
$t\in(0,T]$, where $\varphi_{n}$ and $f_{n}(t)$ are corresponding Fourier
coefficients. Moreover, the series obtained after applying term-wise the
operators $D_{t}^{\rho}$ and $A$ also converge absolutely and uniformly on
$x\in\mathbb{T}^{N}$ and for each $t\in(0,T]$.

######  Remark 1.3.

When $a>\frac{N}{2}$, according to the Sobolev embedding theorem, all
functions in $L^{a}_{2}(\mathbb{T}^{N})$ are $2\pi$-periodic continuous
functions. The fulfillment of the inverse inequality $a\leq\frac{N}{2}$,
admits the existence of unbounded functions in $L^{a}_{2}(\mathbb{T}^{N})$
(see, for example, [9]). Therefore, condition $a>\frac{N}{2}$ for function $f$
of this theorem is not only sufficient for the statement to be hold, but it is
also necessary.

A result similar to Theorem 1.2 was obtained in the recent paper [8] for a
subdiffusion equation with the Riemann-Liouville derivative. But the condition
found for the initial function $\varphi(x)$ in that work is less restrictive.
For example, in the one-dimensional case, it has been proved that it is
sufficient to require $\varphi(x)\in C(\mathbb{T}^{1})$. As will be proved at
the end of this paper, even Holder continuous with exponent $a=\frac{1}{2}$
for the initial function is not sufficient for the validity of Theorem 1.2 in
the one-dimensional case.

In conclusion, note that since our solution is established on the
eigenfunction expansion of elliptic operator, the method proposed in this
article can be used to an arbitrary domain and an elliptic operator with
variable coefficients.

##  2\. Uniqueness

In this section we prove the uniqueness of the solution to problem
(1.1)-(1.2). Note that in order for problem (1.1)-(1.2) to have a unique
solution, it is sufficient that the functions $f$ and $\varphi$ be continuous.

Suppose that initial-boundary value problem (1.1) - (1.2) has two classical
solutions $u_{1}(x,t)$ and $u_{2}(x,t)$. Our aim is to prove that
$u(x,t)=u_{1}(x,t)-u_{2}(x,t)\equiv 0$. Since the problem is linear, then we
have the following homogenous problem for $u(x,t)$:

(2.1) $D_{t}^{\rho}u(x,t)+A\,u(x,t)=0,\quad x\in\mathbb{T}^{N},\quad t>0;$
(2.2) $u(x,0)=0,\quad x\in\mathbb{T}^{N}.$

Let $u(x,t)$ be a solution of problem (2.1)-(2.2). Consider the function

(2.3) $w_{n}(t)=\int\limits_{\mathbb{T}^{N}}u(x,t)e^{inx}dx,\quad
n\in\mathbb{Z}^{N}.$

By virtue of equation (2.1), we can write

$D_{t}^{\rho}w_{n}(t)=\int\limits_{\mathbb{T}^{N}}D_{t}^{\rho}u(x,t)e^{inx}dx=-\int\limits_{\mathbb{T}^{N}}A\,u(x,t)e^{inx}dx,\quad
t>0,$

or, integrating by parts,

$D_{t}^{\rho}w_{k}(t)=-\int\limits_{\mathbb{T}^{N}}u(x,t)A\,e^{inx}dx=-|n|^{2}\int\limits_{\mathbb{T}^{N}}u(x,t)e^{inx}dx=-|n|^{2}w_{n}(t),\quad
t>0.$

Using in (2.3) the homogenous initial condition (2.2), we have the following
Cauchy problem for $w_{n}(t)$:

$D_{t}^{\rho}w_{n}(t)+|n|^{2}w_{n}(t)=0,\quad t>0;\quad w_{n}(0)=0.$

This problem has the unique solution; therefore, the function defined by
(2.3), is identically zero: $w_{n}(t)\equiv 0$ (see for example, [1] p. 17 ,
[10]). From completeness in $L_{2}(\mathbb{T}^{N})$ of the system of
eigenfunctions $\\{e^{inx}\\}$, we have $u(x,t)=0$ for all
$x\in\mathbb{T}^{N}$ and $t>0$. Hence the uniqueness is proved.

##  3\. Existence

Proof of existence based on the following lemma (see M.A. Krasnoselski et al.
[11], p. 453), which is a simple corollary of the Sobolev embedding theorem.

######  Lemma 3.1.

Let $\sigma>1+\frac{N}{4}$. Then for any $|\alpha|\leq 2$ operator
$D^{\alpha}(\hat{A}+1)^{-\sigma}$ (completely) continuously maps from
$L_{2}(\mathbb{T}^{N})$ into $C(\mathbb{T}^{N})$ and moreover the following
estimate holds true

(3.1) $||D^{\alpha}(\hat{A}+1)^{-\sigma}g||_{C(\mathbb{T}^{N})}\leq
C||g||_{L_{2}(\mathbb{T}^{N})}.$

Proof of this lemma see in [8].

One can easily verify that the function (1.4) formally satisfies the
conditions of problem (1.1)-(1.2). In order to prove that function (1.4) is
actually a solution to the problem, it remains to substantiate this formal
statement, i.e. show that the operators $A$ and $D_{t}^{\rho}$ can be applied
term by term to the series (1.4). To do this we remind the following
asymptotic estimate of the Mittag-Leffler function with a sufficiently large
negative argument $t$ and an arbitrary complex number $\mu$ (see, for example,
[12], p. 134)

(3.2) $|E_{\rho,\mu}(-t)|\leq\frac{C}{1+t},\quad t>0.$

We will also use a coarser estimate with a positive $\lambda$ and
$0<\varepsilon<1$:

(3.3) $|t^{\rho-1}E_{\rho,\mu}(-\lambda
t^{\rho})|\leq\frac{Ct^{\rho-1}}{1+\lambda t^{\rho}}\leq
C\lambda^{\varepsilon-1}t^{\varepsilon\rho-1},\quad t>0,$

which is easy to verify. Indeed, let $t^{\rho}\lambda<1$, then
$t<\lambda^{-1/\rho}$ and

$t^{\rho-1}=t^{\rho-\varepsilon\rho}t^{\varepsilon\rho-1}<\lambda^{\varepsilon-1}t^{\varepsilon\rho-1}.$

If $t^{\rho}\lambda\geq 1$, then $\lambda^{-1}\leq t^{\rho}$ and

$\lambda^{-1}t^{-1}=\lambda^{-1+\varepsilon}\lambda^{-\varepsilon}t^{-1}\leq\lambda^{\varepsilon-1}t^{\varepsilon\rho-1}.$

Note the series (1.4) is in fact the sum of two series. Consider the following
partial sums of the first series:

(3.4)
$S^{1}_{k}(x,t)=\sum\limits_{|n|^{2}<k}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,e^{inx},$

and suppose that function $\varphi$ satisfies the condition of Theorem 1.2,
i.e. for some $\tau>\frac{N}{4}$

$\sum\limits_{n\in\mathbb{Z}^{N}}|n|^{4\tau}|\varphi_{n}|^{2}\leq
C_{\varphi}<\infty.$

Since $\hat{A}^{-\tau-1}e^{inx}=|n|^{-2(\tau+1)}e^{inx}$, we may rewrite the
sum (3.4) as

$S^{1}_{k}(x,t)=\hat{A}^{-\tau-1}\sum\limits_{|n|^{2}<k}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,|n|^{2(\tau+1)}\,e^{inx}.$

Therefore by virtue of Lemma 3.1 one has

$||D^{\alpha}S^{1}_{k}||_{C(\mathbb{T}^{N})}=||D^{\alpha}\hat{A}^{-\tau-1}\sum\limits_{|n|^{2}<k}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,|n|^{2(\tau+1)}\,e^{inx}||_{C(\mathbb{T}^{N})}\leq$
(3.5) $\leq
C||\sum\limits_{|n|^{2}<k}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,|n|^{2(\tau+1)}\,e^{inx}||_{L_{2}(\mathbb{T}^{N})}.$

Using the orthonormality of the system $\\{e^{inx}\\}$, we will have

(3.6) $||D^{\alpha}S^{1}_{k}||^{2}_{C(\mathbb{T}^{N})}\leq
C\sum\limits_{|n|^{2}<k}\big{|}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,|n|^{2(\tau+1)}\big{|}^{2}.$

Application of estimate (3.2) and inequality
$|n|^{2}(1+|n|^{2}t^{\rho})^{-1}<t^{-\rho}$ gives

$\sum\limits_{|n|^{2}<k}\big{|}E_{\rho,1}(-|n|^{2}t^{\rho})\varphi_{n}\,|n|^{2(\tau+1)}\big{|}^{2}\leq
Ct^{-2\rho}\sum\limits_{|n|^{2}<k}|n|^{4\tau}|\varphi_{n}|^{2}\leq
Ct^{-2\rho}C_{\varphi}.$

Therefore we can rewrite the estimate (3.7) as

$||D^{\alpha}S^{1}_{k}||^{2}_{C(\mathbb{T}^{N})}\leq Ct^{-2\rho}C_{\varphi}.$

This implies uniformly on $x\in\mathbb{T}^{N}$ convergence of the
differentiated sum (3.4) with respect to the variables $x_{j}$ for each
$t\in(0,T]$. On the other hand, the sum (3.5) converges for any permutation of
its members as well, since these terms are mutually orthogonal. This implies
the absolute convergence of the differentiated sum (3.4) on the same interval
$t\in(0,T]$.

Now we consider the second part of the series (1.4):

(3.7)
$S^{2}_{k}(x,t)=\sum\limits_{|n|^{2}<k}\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})\,d\xi\,e^{inx}$

and suppose that function $f(x,t)$ satisfies all the conditions of Theorem
1.2, i.e. the following series converges uniformly on $t\in[0,T]$ for some
$\tau>\frac{N}{4}$:

$\sum\limits_{n\in\mathbb{Z}^{N}}|n|^{4\tau}|f_{n}(t)|^{2}\leq C_{f}<\infty.$

We choose a small $\varepsilon>0$ in such a way, that
$\tau+1-\varepsilon>1+\frac{N}{4}$. Since
$\hat{A}^{-\tau-1+\varepsilon}e^{inx}=|n|^{-2(\tau+1-\varepsilon)}e^{inx}$, we
may rewrite the sum (3.7) as

$S^{2}_{k}(x,t)=\hat{A}^{-\tau-1+\varepsilon}\sum\limits_{|n|^{2}<k}\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})\,d\xi\,|n|^{2(\tau+1-\varepsilon)}e^{inx}.$

Then by virtue of Lemma 3.1 one has

$||D^{\alpha}S^{2}_{k}||_{C(\mathbb{T}^{N})}=||D^{\alpha}\hat{A}^{-\tau-1+\varepsilon}\sum\limits_{|n|^{2}<k}\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})\,d\xi\,|n|^{2(\tau+1-\varepsilon)}e^{inx}||_{C(\mathbb{T}^{N})}\leq$
(3.8) $\leq
C\big{|}\big{|}\sum\limits_{|n|^{2}<k}\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})\,d\xi\,|n|^{2(\tau+1-\varepsilon)}e^{inx}\big{|}\big{|}_{L_{2}(\mathbb{T}^{N})}.$

Using the orthonormality of the system $\\{e^{inx}\\}$, we will have

$||D^{\alpha}S^{2}_{k}||^{2}_{C(\mathbb{T}^{N})}\leq
C\sum\limits_{|n|^{2}<k}\big{|}\int\limits_{0}^{t}f_{n}(t-\xi)\xi^{\rho-1}E_{\rho,\rho}(-|n|^{2}\xi^{\rho})\,d\xi\,|n|^{2(\tau+1-\varepsilon)}\big{|}^{2}.$

Now we use estimate (3.3) and apply the generalized Minkowski inequality. Then

$||D^{\alpha}S^{2}_{k}||^{2}_{C(\mathbb{T}^{N})}\leq
C\bigg{(}\int\limits_{0}^{t}\xi^{\varepsilon\rho-1}\big{(}\sum\limits_{|n|^{2}<k}|n|^{4\tau}|f_{n}(t-\xi)|^{2}\big{)}^{1/2}d\xi\bigg{)}^{2}\leq
C\cdot C_{f},$

where $C$ depends on $T$ and $\varepsilon$. Hence, using the same argument as
above, we see that the differentiated sum (3.7) with respect to the variables
$x_{j}$ converges absolutely and uniformly on
$(x,t)\in\mathbb{T}^{N}\times[0,T]$.

Further, from equation (1.1) one has

$D_{t}^{\rho}\big{(}S^{1}_{k}+S^{2}_{k}\big{)}=-A\big{(}S^{1}_{k}+S^{2}_{k}\big{)}+\sum\limits_{|n|^{2}<k}f_{n}(t)e^{inx}.$

Absolutely and uniformly convergence of the latter series can be proved as
above.

Thus Theorem 1.2 is completely proved.

##  4\. Counterexample

In this section, we will discuss the importance of the condition
$a>\frac{N}{2}$ of Theorem 1.2. When $a>\frac{N}{2}$ all functions in
$L_{2}^{a}(\mathbb{T}^{N})$ belong to $C(\mathbb{T}^{N})$. As noted in Remark
1.3, if $a=\frac{N}{2}$, then in the class $L_{2}^{a}(\mathbb{T}^{N})$ there
exist unbounded functions, as a consequence of which problem (1.1)-(1.2)
certainly does not have classical solutions. Therefore, the question naturally
arises: is it possible to replace, for example, condition $\varphi\in
L_{2}^{a}(\mathbb{T}^{N})$, $a>\frac{N}{2}$, of Theorem 1.2 by condition

(4.1) $\varphi\in L_{2}^{\frac{N}{2}}(\mathbb{T}^{N})\cap C(\mathbb{T}^{N})?$

The following example answers this question in the negative.

Let $N=1$. In the class of periodic functions, we seek a solution to the
following problem

(4.2) $\left\\{\begin{aligned} &D_{t}^{\rho}u(x,t)-u_{xx}(x,t)=0,\quad 0<t\leq
T;\\\ &u(x,o)=\varphi(x).\end{aligned}\right.$

If $\varphi$ satisfies the condition of Theorem 1.2, then the unique solution
of the problem has the form

(4.3) $u(x,t)=\lim\limits_{k\rightarrow\infty}\sum\limits_{|n|\leq
k}\varphi_{n}E_{\rho,1}(-|n|^{2}t^{\rho})\,e^{inx}.$

Recall, that the classes $C^{a}(\mathbb{T}^{1})$ are usually defined as
follows: $2\pi$-periodic function $\varphi\in C^{a}(\mathbb{T}^{1})$ if and
only if

$|\varphi(x)-\varphi(y)|\leq C|x-y|^{a}.$

Consider the following function, first studied by Hardy and Littlewood (see
[13], proof of Theorem (3.10)):

$\Phi(x)=\sum\limits_{n=1}^{\infty}\frac{e^{in\ln n}}{n}e^{inx}.$

The real and imaginary parts of this function belong to the class
$C^{\frac{1}{2}}(\mathbb{T}^{1})$ (see [13], Chapter V, paragraph 4). Set
$\varphi(x)=\Re(\Phi(x))$, where $\Re(z)$ is a real part of the complex number
$z$. Then $\varphi\in C^{\frac{1}{2}}(\mathbb{T}^{1})$ and it is not hard to
see, that

$\sum\limits_{n=1}^{\infty}\sqrt{(\varphi^{c}_{n})^{2}+(\varphi^{s}_{n})^{2}}=+\infty,$

where $\varphi^{c}_{n}$ are the coefficients of the function $\varphi(x)$ in
$\cos nx$ and $\varphi^{s}_{n}$ \- in terms of $\sin nx$. Obviously, the
function $\varphi(x)$ also belongs to the class (4.1) and if we denote
$\varphi_{n}=\frac{1}{2}(\varphi^{c}_{n}-i\varphi^{s}_{n})$ and
$\varphi_{-n}=\overline{\varphi_{n}}$, then

(4.4) $\sum\limits_{n\in\mathbb{Z}^{1}}|\varphi_{n}|=+\infty.$

Suppose that the solution to problem (4.2) has the form (4.3). Let us show
that the series (4.3) differentiated twice with respect to the variable $x$
does not converge absolutely, that is, for the function $\varphi(x)$, defined
above, the statement of Theorem 1.2 does not hold. Indeed, set

$(u_{k_{0}})_{xx}(x,t)=-\sum\limits_{k_{0}\leq|n|}\varphi_{n}E_{\rho,1}(-|n|^{2}t^{\rho})\,|n|^{2}\,e^{inx}.$

In order for this series to converge uniformly and absolutely with respect to
$x\in\mathbb{T}^{1}$ and for each $t\in(0,T]$ it is necessary that the number
series

$U_{k_{0}}=\sum\limits_{k_{0}\leq|n|}|\varphi_{n}E_{\rho,1}(-|n|^{2}t^{\rho})|\,|n|^{2}$

converge for some $t>t_{0}>0$. Now we remind the following asymptotic estimate
of the Mittag-Leffler function with a sufficiently large negative argument

$E_{\rho,1}(-|n|^{2}t^{\rho})=\frac{1}{\Gamma(1-\rho)}\cdot\frac{1}{|n|^{2}t^{\rho}}+O\bigg{(}\frac{1}{|n|^{2}t^{\rho}}\bigg{)}^{2}.$

Hence, for sufficiently large $k_{0}$, we have

$U_{k_{0}}=\frac{1}{t^{\rho}\Gamma(1-\rho)}\cdot\sum\limits_{k_{0}\leq|n|}|\varphi_{n}|+O(1),$

where $O(1)$ depends on $t_{0}$ and $\rho$. Since (4.4) this series does not
converge.

Thus condition $\varphi\in L_{2}^{a}(\mathbb{T}^{N})$, $a>\frac{N}{2}$, of
Theorem 1.2 cannot be replaced by the condition (4.1) in at least one-
dimensional case.

In conclusion, note that a similar result with Theorem 1.2 for equation (1.1)
with the Riemann-Liouville derivative is valid for all functions $\varphi\in
C(\mathbb{T}^{N})\cap L_{2}^{a-2}(\mathbb{T}^{N})$ (see [8]).

## References

  * [1] Pskhu, A.V.: Fractional partial differential equations (in Russian), M. NAUKA (2005)
  * [2] _Kochubei A., Luchko Yu._. Handbook of Fractional Calculus with Applications. Volume 2: Fractional Differential Equations. De Gruyter 2019.
  * [3] Kilbas, A.A., Srivastava, H.M., Trujillo, J.J.: Theory and applications of fractional differential tquations, ELSEVIER (2006)
  * [4] Yu. Luchko, Some uniqueness and existence results for the initial-boundary value problems for the generalized time-fractional diffusion equation, Comput. Math. Appl. 59 (2010) 1766-1772.
  * [5] Ruzhansky, M., Tokmagambetov, N., Torebek, B.T.: On non-local problem for a multi-term fractional diffusion-wave equation, arXiv:1812.01336v2 [math. AP] 5 Dec (2018)
  * [6] K. Sakamoto and M. Yamamoto, Initial value/boundary value problems for fractional diffusion-wave equations and applications to some inverse problems, J. Math. Anal. Appl. 382 (2011) 426-447.
  * [7] R. Ashurov, O. Muhiddinova, Initial-boundary value problem for a time-fractional subdiffusion equation with an arbitrary elliptic differential operator, Lobachevski Journal of Mathematics, 2021, Vol. 42, No. 3, pp. 517-525.
  * [8] R. Ashurov, O. Muhiddinova, Initial-boundary value problem for a time-fractional subdiffusion equation on the torus,https://arxiv.org/abs/2105.07415
  * [9] Alimov, Sh.A., Ashurov, R.R., Pulatov, A.K.: Multiple Fourier Series and Fourier Integrals. Commutative Harmonic Analysis, Springer, Berlin (1992)
  * [10] R. Ashurov, A. Cabada and B. Turmetov, “Operator method for construction of solutions of linear fractional differential equations with constant coefficients”, Fractional calculas and applied analysis, 1, 229-252 (2016).
  * [11] Krasnoselski, M.A., Zabreyko, P.P., Pustilnik E.I., Sobolevski, P.S.: Integral operators in the spaces of integrable functions (in Russian), M. NAUKA (1966)
  * [12] M. M. Dzherbashian [=Djrbashian], _Integral Transforms and Representation of Functions in the Complex Domain_ (in Russian), M. NAUKA, 1966.
  * [13] A. Zygmund, _Trigonometric series_ , Vol. 2, Cambridge, The univereity Press (1959).