{'fraction_non_alphanumeric': 0.04154125909204916, 'fraction_numerical': 0.016114873338349635, 'mean_word_length': 4.548269264219864, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Neural Network Approach for Stochastic Optimal Control

We present a neural network approach for approximating the value function of high-dimensional stochastic control problems. Our training process simultaneously updates our value function estimate and identifies the part of the state space likely to be visited by optimal trajectories. Our approach leverages insights from optimal control theory and the fundamental relation between semi-linear parabolic partial differential equations and forward-backward stochastic differential equations. To focus the sampling on relevant states during neural network training, we use the stochastic Pontryagin maximum principle (PMP) to obtain the optimal controls for the current value function estimate. By design, our approach coincides with the method of characteristics for the non-viscous Hamilton-Jacobi-Bellman equation as the dynamics become deterministic. Our training loss consists of a weighted sum of the objective functional of the control problem and penalty terms that enforce the HJB equations along the sampled trajectories. Importantly, training is unsupervised in that it does not require solutions of the control problem. In our numerical experiments, we compare our method to existing solvers for a more general class of semi-linear PDEs. Using a two-dimensional model problem, we demonstrate the importance of the stochastic PMP to inform the sampling. For a 100-dimensional benchmark problem, we demonstrate that our approach improves accuracy and time-to-solution.

1. Introduction. We consider Stochastic Optimal Control (SOC) problems that require finding a policy to control randomly perturbed dynamical systems so as to optimize a given objective functional. Problems of this type arise in many fields including finance, biology, robotics and many other engineering applications; see, e.g., [10,34] for references and extensive theoretical discussion on SOC problems.
Dynamic programming (DP) is a prominent framework for solving SOC problems. At its core, DP seeks to find the value function, which assigns every state of the system the optimal cost-to-go. For many problems, recovering the optimal control from the value function is straightforward. One way to compute the value function is by solving the Hamilton-Jacobi-Bellman (HJB) equation [2,37], which is a semi-linear parabolic Partial Differential Equation (PDE). The two main challenges in solving the HJB equations is the forward-backward structure and the high dimensionality.
The space dimension of the HJB equation equals the state dimension of the dynamical system to be controlled. Hence, many numerical schemes for solving PDEs cannot be applied to realistic problem instances due to the Curse-of-Dimensionality (CoD); for example, the computational costs of approaches that discretize the state space with a mesh typically grows exponentially in the dimensions of the system. Hence, the applicability of existing approaches that employ spatial discretizations (for example, [4,6,16,20,32,22,23,32,35]) is limited to state dimensions d ≤ 3. Mesh-free methods based on radial basis function methods (for example, [36]) can be effective for slightly larger d but ultimately also suffer from CoD for d = O(100), as we consider here.
A first requirement for mitigating the CoD in SOC problems is to alleviate the need for a spatial discretization. This can be achieved by using a nonlinear version of the Feynman-Kac lemma and replacing the HJB equation by a system of Forward-Backward Stochastic Differential Equations (FBSDEs); see, for example, [5,28,37]. Discretizing the FBSDE system using the classical Euler Maruyama scheme [19] yields mesh-free numerical schemes that have gained importance for SOC problems and nonlinear second-order PDEs; see, for example, [3,8,9,13,25] and references therein.
A second requirement for mitigating the CoD is to parameterize the value function effectively in high-dimensions. Due to their universal approximation property and many advances in deep learning, the use of neural networks (NNs) as function approximators has been gathering significant attention recently. The idea of combining the BSDE approach with NNs has been pioneered by the seminal works [7,12] that use neural networks to solve high-dimensional semi-linear parabolic PDEs and include an example for HJB equations as a special case. Around the same time, [33] achieves similarly promising results for neural-network based solution of the HJB equation. These works developed learning algorithms for approximating the value function or its gradient using the FBSDE system. Their success has sparked several extensions that improve the NN architectures and consider other types of PDEs; see, for example, [15,17,29].
Even after combining FBSDEs and deep learning, solving the HJB equation globally remains affected by the curse of dimensionality as it would require sampling the entire state-space. Fortunately, the optimal state trajectories of many SOC problems do not cover the entire space. This motivates us to develop semi-global approaches, that is, we seek to reliably estimate the value function in those parts of the state space that are likely to be visited when following an optimal policy. Clearly, these tasks are interrelated: Finding the relevant parts of the state space depends on the value function and, vice versa, the value function needs to be trained using samples from the state space.
The key idea of our approach is to use the stochastic Pontryagin maximum principle (PMP) [31] to define the forward system in the FBSDE approach in terms of the value function. The PMP provides a set of necessary optimality conditions, known as the Hamiltonian system, and a feedback form that allows recovering the optimal control from the value function; see, for example, [37]. As we discuss in more detail later, the choice of the forward SDE is crucial as it is used to sample the parts of the state space along which the value function estimate is improved using a loss that includes the backward SDE. It is worth noting that the FBSDE framework allows choosing the forward system almost arbitrarily. Thus, our choice is theoretically on par with the standard Brownian motion used for exploring the state space in [12,33], but as our numerical experiments suggest enable solving a wider class of practical problems; see section 4.
A theoretical advantage of our proposed scheme over existing approaches is its consistency with deterministic optimal control problems. In the absence of uncertainty in the dynamics, our approach reduces to the method of characteristics for the HJB equation. Hence our approach can also be seen as an extension of the neural network approaches for deterministic optimal control proposed in [27,21].
The rest of the paper is organised as follows: In section 2 we introduce the parts of SOC theory used to obtain our proposed forward SDE. In section 3, we provide the NN approach to learn the value function using the FBSDE reformulation from section 2. Using an intuitive 2D example, we show the importance of sampling. To illustrate the potential of our method, we consider the 100-dimensional benchmark problem also used in [12,33], propose a modification that leads to a more challenging sampling problem, and investigate the scaling of our approach with respect to problem dimension. Finally we conclude the paper and discuss future directions.
2. Stochastic Optimal Control Background. In this section, we describe the Stochastic Optimal Control (SOC) problems considered in this work and review the key results from SOC theory that motivates our approach; our discussion follows [37] and we refer to this text book for a more comprehensive background. We first introduce the SOC problem and then review its underlying theory; to be specific, we review the stochastic Pontryagin Maximum Principle (PMP), the Hamilton-Jacobi-Bellman (HJB) equation, and its reformulation into a system of forward backward stochastic differential equations (FBSDEs) obtained from a nonlinear version of the Feynman-Kac formula.

Stochastic Optimal Control
Problem. Let (Ω, F, F = {F t } t≥0 , P) be a given complete probability space, W (s) be a d-dimensional Brownian motion on (Ω, F, F, P) where s denotes the time. For a fixed initial state x at some time 0 < t < T < ∞, we seek to control the randomly perturbed dynamical system Here, z t,x : [t, T ] → R d describes the state and u t,x : [t, T ] × R d → U describes the control of the system, the function σ : [t, T ] × R d → R d×d represents the diffusion coefficient, and f : [t, T ] × R d × U → R d represents the drift of the system. We assume that the set of admissible controls U ⊂ R k is closed. We seek to minimize the objective functional which is comprised of the running cost L : [t, T ] × R d × U → R and the terminal cost G : R d → R. Here, the expectation is taken with respect to perturbances of the dynamics (2.1) that is described by the probability measure P. We assume sufficient regularity conditions on f , σ, G and L, see [37, Chapter 2] for a list of assumptions. The value function assigns the optimal cost-to-go to any initial state, that is, and a solution u * t,x to (2.3) incurring this minimum value is called an optimal control. The generalized Hamiltonian, H : To make it notationally convenient, in the rest of the paper, we drop the second argument for the controls and denote controls by u t,x (s).

Stochastic Pontryagin Maximum Principle.
The PMP provides firstorder necessary conditions for the SOC problem and also states that the optimal control u * t,x must satisfy an (extended) Hamiltonian system along the optimal state and adjoint trajectory. This is made precise by the following result from [ where s ∈ [t, T ] and the optimal control satisfies for almost all s ∈ [t, T ], P-almost surely.
We note that the optimal control defined in (2.6) only depends on the adjoint variable p t,x but not on M t,x since the diffusion σ(·, ·) does not depend on the control in our case.
We further assume that there exists a unique continuous closed form solution to (2.6). Although not demonstrated in this work, this assumption can be weakened to include implicitly defined functions as long as they can be obtained efficiently; this allows, for example, modeling more general convex running costs.
We note that when the control satisfies (2.6), the dynamics in (2.1) is equal to (2.7) The system of equations (2.5)-(2.7) is called the stochastic Hamiltonian system, where the maximum condition (2.6) corresponds to the variational inequality for the control.
Finding a tuple (z * t,x , u * t,x , p t,x , M t,x ) that satisfies the PMP can be extremely difficult. However, when the value function Φ is differentiable, (p t,x , M t,x ) satisfying (2.5) can be obtained from Φ; this is formalized in the following theorem that, with weaker assumptions, can be found in [37,Chapter 5].
Theorem 2.2 simplifies (2.6) to This relation along with (2.7) is one of the key ingredients of our numerical solution approach. Equation (2.9) characterizes optimal control in a feedback or closed-loop form, which is of utmost importance in many real-life applications. This is because optimal control can be quickly computed at any given point in time and space, when the value function is known and its gradient is readily available. Thereby, avoiding re-computation of optimal controls for cases when multiple evaluations are needed for different time or states.

Hamilton-Jacobi-Bellman Equation.
To help approximate the value function Φ, we also use the fact that Φ satisfies the Hamilton-Jacobi-Bellman (HJB) equation, which is a result of the Dynamic Programming (DP) method or Bellman's principle. We state the following result taken from [37] (2.10) The smoothness of Φ can be relaxed to continuity in the weaker sense of viscosity solutions [37,Section 5,Chapter 4].
Using the relation between PMP and DP from Theorem 2.2 the HJB equation, (2.10), can now be re-written as the following second order parabolic PDE: In addition, by the envelope theorem, it follows that ∇ p H = ∇ p H and ∇ M H = ∇ M H. This simplifies the computation of optimal trajectories, which can now be expressed via the value function as: (2.12) These modified dynamics do not explicitly involve the control, which reduces the problem solely to state variable. This idea was proposed for deterministic optimal control problems in [27].

FBSDE Formulation.
One way to avoid the need for a spatial discretization of the HJB equation (2.11) is to use a non-linear version of Feynman-Kac formula and obtain an equivalent system of stochastic differential equations. This idea has been applied to a variety of nonlinear parabolic/elliptic PDEs; see, e.g., [1,24,5,28,30].
Our FBSDE system uses (2.12) as the forward system to sample trajectories. Along those trajectories, we note that the solution to the HJB equation (2.11) must satisfy the backward SDE dΦ(s, z(s)) = ∇ z Φ(s, z(s)) σ(s, z(s)) dW (s) − L(s, z(s), u * (s))ds, Φ(T, z(T )) = G(z(T )). (2.13) It is important to stress that our choice of the forward system (2.12) is the key difference to other existing works. For example, [12] and [33] use the standard Brownian motion. While both of these choices lead to a valid FBSDE system for (2.11), we advocate for including the control in the dynamics as motivated by stochastic PMP (2.7). As our numerical experiments demonstrate, focusing the sampling along optimal trajectories can lead to more accurate and efficient value function approximations.
3. Neural Network Approach. In this section, we present a neural network framework for approximating the value function of the stochastic optimal control problem defined by the objective functional (2.2) and dynamics (2.1). The theoretical foundation of our framework are given by the PMP, FBSDE, and Dynamic Programming as presented in the previous section. The key idea is to approximate the value function Φ in (2.3) by a neural network and compute the control using the feedback form (2.9). What sets our framework apart is the use of the feedback form to guide the sampling during training. Thereby we seek to learn to explore the relevant part of the state space. We also derive and experiment with various loss functions that are based on the control objective (2.2), the BSDE (2.13), and the HJB (2.11).
3.1. Neural Network Approximation. The first building block of our framework is to parameterize the value function using a neural network. Since finding an effective network architecture for any learning task is both crucial and an open research topic, we treat this as a modular component. Our framework can be used with any scalar-valued neural network that takes inputs in R d+1 as long as it is twice differentiable with respect to its inputs.
Among the networks we use in our numerical experiments is the Multi-Layer Perceptron (MLP) model used in [33]. As an alternative, which also satisfies the regularity needed, we propose the residual network also used for deterministic control in [27]. The network is given by θ N N , A, b, c). Here the inputs y = (s, z(s)) ∈ R d+1 correspond to time-space, N N (y; θ N N ) : R d+1 → R m is a neural network, and θ contains the trainable weights: where rank γ= min(10, d + 1) limits the number of parameters in A A. Here, A, b, and c model quadratic potentials, i.e., linear dynamics; N N models nonlinear dynamics.
In our experiments, for N N , we either use a MLP [11] a 0 = act(K 0 y + b 0 ) or a residual neural network (ResNet) [14] a 0 = act(K 0 y + b 0 ) with neural network weights θ N N =(K 0 , . . . , , and {K 1 , . . . , K M } ∈ R m×m with M being the depth of the network. The choice of the element-wise nonlinearity act(·) is discussed in the respective experiments.
3.2. Training Problem. Ideally, we would choose θ such that Φ(s, z; θ) is equal to the value function of the control problem globally, that is, for all (s, z) ∈ [t, T ]×R d . Since this is known to be cursed by the dimensionality for reasonable problem sizes, we resort to a semi-global approach which enforces this property at randomly sampled points in the space time domain.
To generate samples, we first obtain initial states x ∼ ρ from some (possibly Dirac) distribution ρ and then use an Euler Maruyama scheme with N + 1 equidistant time points s 0 , . . . , s N and step size ds = (T − t)/N . This yields a state trajectory starting at z 0 = x via is the optimal control obtained from the feedback, that is, form (2.6) A few comments are in place. First, it is important to note that due to the feedback form, the sampled trajectories depend on the parameters of the value function.
Second, the addition of this drift term, motivated by control theory, is the key difference to neural network solvers for the more general class of semi-linear elliptic PDEs [12,33]. Third, the drift term can also be motivated by the fact that for σ → 0, the trajectories above approximate the characteristic curves of the non-viscous HJB equation; thereby our SOC approach coincides with that for deterministic OC in [27]. To simplify our notation below, we collect the states, control, and noise along the discrete trajectories in (3.4) column-wise in the matrices To learn the parameters of the neural networks in an unsupervised way (that is, assuming neither analytic values of Φ nor optimal control trajectories), we approximately solve the minimization problem where the terms in the objective function consist of penalty functions for violations of the BSDE system and the HJB PDE, the control objective, and penalty terms for the terminal condition, respectively, and are defined below. The relative influence of each term is controlled by the components of β ∈ R 5 + . Different choices of β allow us to experiment with different learning approaches; for example, setting β 1 = β 4 = β 5 = 1 and β 2 = β 3 = 0 provides the same loss function as in [33] while β 1 = 0 and β i > 0, i ∈ {2, 3, 4, 5} gives the loss function used for deterministic OC problems in [27]. We penalize the violation of the BSDE (2.13) via (3.6) where we use the abbreviations Φ i (θ) := Φ(s i , z i ; θ) and ∇Φ i (θ) := ∇Φ(s i , z i ; θ) and the exponent p ∈ {1, 2}. Similarly, the HJB penalty term reads . Finally, we approximate the objective functional via In principle, any stochastic approximation approach can be used to approximately solve the above optimization problem. Here, we use Adam and sample a minibatch of trajectories originating in n i.i.d. samples from ρ.

Implementation.
We implement and test our proposed approach in two software environments.
To obtain a direct comparison with [33] we modify the FBSNN code accompanying the paper. To this end, we created a publicly available fork at https://github.com/ EmoryMLIP/FBSNNs. Our main modifications are adding the control objective in the training loss and also adding the proposed drift to the forward dynamics. The experiments carried out here use a standard MLP as the neural network without the quadratic term added in (3.1).
In order to further simplify the experimentation, we also implement our own PyTorch code available at https://github.com/EmoryMLIP/NeuralHJB. Our implementation contains all loss terms in (3.5). We implement both sampling techniques: pure random walk and the proposed one informed by PMP. This facilitates comparisons of our approach with other available methods and simplifies developing new examples.
We tested most of our examples using either Intel Xeon E5-4627 CPU or Nvidia P100 GPU.

Numerical Experiments.
We test the efficacy of our proposed algorithm on two different SOC problems. In subsection 4.1, we introduce a two-dimensional trajectory planning problem to visualize the difference between purely random exploration and our proposed sampling scheme. In subsection 4.2, we compare our approach to those in [7,12] using a 100-dimensional benchmark problem. For the original version of this problem, our method shows faster initial convergence and time-to-solution but comparable accuracy. We modify the terminal cost of this problem to further highlight the importance of the feedback form in the sampling (see subsection 4.2.2) and experimentally show that the growth of network size is not exponential with increasing dimensions.

2D Trajectory Planning Problem.
To visualize the behaviour of our PMP-based sampling approach, we consider a two-dimensional test problem.
The problem consists of planning an optimal trajectory from the initial state x ∼ ρ = N ((−1.5, −1.5) , I 2 ) to the target x target = (1.5, 1.5) . To make the problem interesting, a hill is placed at the origin, which adds height-dependent cost for traveling around that region given by The dynamics for the problem read The running cost and terminal cost of the problem are given respectively, by The corresponding HJB equation is With the above calculations, the sampling SDE (3.4) simplifies to We approximate the value function using the model in (3.1) where N N is a twolayer residual neural network with 32 neurons per layer and the rank of the quadratic term is γ = 3. Overall, the model has 1229 trainable weights. As activation function we chose the anti-derivative of the tanh function as in [27]. We use the penalty parameters β = (1, 1, 1, 0.1, 0.1), that is, we enable both the penalty terms, P BSDE and P HJB . To approximately solve (3.5) we use 3000 steps of the Adam optimizer and batch size of 64. After 1500 iterations, we reduce the learning rate from 0.01 to 0.001. The average cost per iteration is about 0.1s using the CPU.
In Figure 1, we visualize the estimate of the value function and the optimal control policy at t = 0 as well as example trajectories for randomly chosen initial states after training the model. As expected, the trajectories travel from the initial points to the target while avoiding the obstacle in the center of the domain.
To gain more insight into the sampling, we store all states visited during training and plot them as two-dimensional histograms for different time points (left to right) in Figure 2. We compare the proposed PMP-based sampling (Figure 2b) to the pure random walk (i.e., without drift) used in [12,33] (Figure 2a). As to be expected, sampling without drift leads to visiting points around the initial states in all (even suboptimal) directions. It can also be seen that there are almost no samples close to the target. With the drift term and after training, the sampled states cover more likely paths between the initial and target states.
Another way to interpret the histogram plots in Figure 2 is that they show the semi-global nature of neural network approaches for SOC problems and, crucially, that the regions of the state space visited in training can vary a lot. Since the loss function in (3.5) penalizes the HJB and BSDE loss only around points sampled from the forward SDE, one would expect that the trained model is more reliable in regions that are frequently visited. (a) Training samples with pure random walk as FSDE as also used in [12] and [33].  At six time points (left to right), we visualize the sampled states as two-dimensional histrograms. As expected, the pure random walk explores the area around the initial state in all (even suboptimal) directions, while the proposed approach learns to sample around approximately optimal trajectories. 4.2. 100-dimensional example. We consider the 100-dimensional benchmark SOC problem also used in [7,12] with initial state x = (0, 0, . . . , 0) ∈ R 100 corresponding to time t = 0. The drift and diffusion of the system are given by Using the first-order necessary condition we get and using this closed form for u, the Hamiltonian is given by Hence, the HJB equation satisfied by the value function, Φ(·, ·), reads and its solution is given by which we use to test the performance of our method. Finally, we note that the forward SDE (3.4) we propose to use for sampling the state space simplifies to

The importance of sampling.
To demonstrate the impact of using the feedback form to sample the state space, we use the same neural network model as in [33], which is given by a 5-layer feed forward neural network with 256 neurons per hidden layer to approximate the solution Φ(s, z). We partition the time interval [0, 1] using 50-uniformly spaced points. We use the same penalty parameters as in the original code, that is, β = (1, 0, 20, 1, 1). We use the Adam optimizer [18] to update the parameters of the network with a batch size of 64 using 50k iterations. The average cost per 100 iterations was 27s using the CPU. For the following experiments, we use (3.5) excluding P HJB penalty and compare our method with FBSNNs in [33].
Ours 0.48% 0.0083% 0.39% 0.012% Table 1: Relative errors for (4.6) obtained using our method and method in [33] In Figure 3 we plot the exact solution (black-dashed line) (4.7), the learned solution using our approach (blue-solid line) and the solution learned using FBSNNs  Fig. 3: Solution to (4.6) obtained using our method and method in [33] (red-solid line) along five random trajectories. In the top row, we present the results obtained after training the networks for 20K iterations with a learning rate of 10 −3 and the bottom row presents the results after training the networks for 20K and 30K iterations with learning rates 10 −3 and 10 −4 , respectively. These results suggest that our approach is approximating the value function better, especially in early iterations, as compared to the FBSNNs.
In Table 1, we also compare the learned solutions to the exact solution Φ in (4.7) by computing the average relative errors, for ten random trajectories. Our method attains lower errors, especially for the initial values and at the earlier iterations.

Shifted target.
In the example above, the minimizer of the terminal function coincides with the initial state x = (0, 0, . . . , 0) . Therefore, even a random walk without drift (as used in [33,12]) will sample around the optimal terminal state, which is critical to accurately approximate the value function. This also means that after training using our approach, the drift term in the sampler is relatively small and that the above experiment does not fully show the advantages of our method.
To shed more light on the importance of sampling, we modify the terminal cost to with z target = [3, 3, . . . , 3] T , so that the target for the state variable z at final time T no longer coincides with the initial state. Similar to the two-dimensional test problem in subsection 4.1, solving the modified problem now requires sampling around the target and we expect to benefit from the added drift term. We compare our method to FBSNNs on the modified problem keeping the same network structure and hyper-parameters. We use a smaller σ = 2 √ 2 5 to improve training speed. We evaluate the performance of the methods using the objective functional J defined in (2.2) at the control obtained from the feedback form via the respective value function approximations. For this experiment we use GPU to train and the results of this comparison are shown in Figure 4.
To reduce the effect of the Brownian motion, we run the experiments for each method on the same problem five times and plot the average values corresponding to training iterations, we also further reduced the effect of the random walk in evaluation for clearer visualization.
As can be seen in Figure 4 (left), our method not only yields faster initial convergence but also achieves a considerably lower control objective. This indicates that the controls obtained from our approach are more effective, that is, they are closer to optimal. Since FBSNNs use a Brownian motion with no drift, the sampling is unlikely to discover the target. Consequently, the generated trajectories in Figure 4 (right) from our method approximately reach the target, while the trajectories obtained from the FBSNN method stay closer to the initial state.

Scalability with problem dimension.
As some evidence that the proposed NN approach can mitigate the curse of dimensionality, we vary the dimension in {20, 40, 60, 80, 100} in subsection 4.2.2. We use the same network architecture from the previous session. For each dimension, we increase the width of the network until the training yields a network whose trajectories end sufficiently close to the target.
With our method we notice that the number of weights required to solve the problem does not increase exponentially with dimensionality, suggesting that our method does ease the CoD. It is however worth mentioning that this experiment does not aim to determine the minimum number of weights required to solve each problem.
5. Discussion. We propose a neural network approach for approximately solving Hamilton-Jacobi-Bellman PDEs arising in high-dimensional stochastic optimal control. Similar to existing approaches [33,12], we parameterize the value function with a neural network and we experiment with different losses to train the weights of the network. What sets our work apart from existing approaches is the use of feedback form given by the stochastic Pontryagin maximum principle to design the forward SDE used to explore the state space during training.
Using an intuitive two-dimensional test problem, we visualize that the improved sampling strategy allows us to effectively learn the value function and determine the relevant regions of the state space; see subsection 4.1. Based on this insight, we modify the 100-dimensional test problem also used in [33,12] by shifting the minimizer of the terminal costs; see subsection 4.2.2. Thereby, we demonstrate that our proposed method dramatically improves the quality of the obtained control.
In additional numerical experiments, we show that the modified forward SDE and the control objective can lead to faster initial convergence compared to the approaches in [33,12] (see subsection 4.2) and show that, for the problem at hand, the network size needed to achieve a given tolerance does not need to grow exponentially.
We refer to our method as a semi-global method for solving the HJB equation since we do not aim at estimating the value function accurately globally. Instead, we seek to approximate the value function well for states likely to be visited by optimal trajectories of the SOC problem. In theory, any point in the state space has a positive probability to be visited due to the stochasticity in the dynamics. However, the histogram plots in Figure 2 suggest that the density of the optimal trajectory is concentrated in a small subset of the state space. Therefore, focussing the exploration to this subsets may have practical advantages over FSDEs that are pure random walks.
Another benefit of our proposed forward SDE compared to purely random exploration is that it coincides with the characteristic curves of the HJB equation as the stochasticity of the system is reduced. Therefore, our work can be seen as an extension of the neural network approaches for deterministic control problems in [27].
Compared to neural network approaches for semi-linear elliptic/parabolic PDEs such as [33,12] it is important to highlight that our approach is only applicable to HJB equations arising in stochastic optimal control. Since our forward SDE is derived from optimality principles, extending it to other high-dimensional PDEs (e.g., Black Scholes and Allen Cahan equations) is not obvious and may be impossible.
In future work, we will apply our approach to problems with non-constant and non-scalar diffusion coefficients. Although the theoretical framework supports this, we are not aware of any practical algorithms for this case. To overcome the challenge of efficiently computing Hessians of the value function needed for HJB penalties, we will use the hessQuik package [26].