Dynamic Programming for Finite Ensembles of Nanomagnetic Particles

We control a ferromagnetic N-spin system which is exposed to thermal fluctuations via an exterior forcing \(\mathbf{u} := (u_1, \ldots , u_N): [0,T] \times \Omega \rightarrow ({{\mathbb {R}}}^3)^N\). A relevant application includes data storage devices, for which it is crucial to control the dynamics of the magnetization \(\mathbf{m} := (m_1, \ldots , m_N): [0,T] \times \Omega \rightarrow ({\mathbb S}^2)^N\) in order to ensure a reliable transport of data which are represented by those magnetic structures. Besides the external force \(H_\mathtt{ext,i}(u_i) = C_\mathtt{ext}\,u_i\), the i-th spin of the ensemble with magnetization \(m_i\) is exposed to the different forces \(H_\mathtt{ani,i}(m_i)\), \(H_\mathtt{d,i}(m_i)\), and \(H_\mathtt{exch,i}(\mathbf{m})\):

• the anisotropic force \(H_\mathtt{ani,i}(m_i) = \mathtt{A} m_i\) for some \( \mathtt{A}\in {{\mathbb {R}}}^{3 \times 3}_\mathtt{diag}\), which favors alignment of \(m_i\) with the given material dependent easy axis \(e \in {{\mathbb {R}}}^3\),

• the ‘stray-field force’ \(H_\mathtt{d,i}(m_i) = - { \mathtt B}_i m_i\) for some \(\mathtt{B}_i \in {{\mathbb {R}}}^{3 \times 3}_\mathtt{diag}\),

• the exchange force \(H_\mathtt{exch,i} (\mathbf{m})\), which penalizes non-alignment of neighboring magnetizations via \(H_\mathtt{exch,i}(\mathbf{m})= -(\mathbf{J}{} \mathbf{m})_i\), for some positive semi-definite \(\mathbf{J} \in {{\mathbb {R}}}^{3N \times 3N}_\mathtt{sym}\).

For every \(1 \le i \le N\), we denote their superposition

as the effective field. The dynamics of \(\mathbf{m}\) for times [0, T] is then governed by the following SDE system (\(1 \le i \le N\)):

Here, \(\mathbf{W} = (W_1, \ldots ,W_N)\) is a \(({{\mathbb {R}}}^3)^N\)-valued Wiener process on a filtered probability space \((\Omega , {\mathcal F}, \{ {{\mathcal {F}}}_t \}_t, {{\mathbb {P}}})\) satisfying the usual conditions to represent thermal fluctuations from the surrounding heat bath. The leading term in the drift part in (1.2) causes a precessional motion of \(m_i\) around \(H_\mathtt{eff,i} (\mathbf{m}, \mathbf{u})\), while the dissipative second term scaled by \(\alpha >0\) favors a time-asymptotic alignment of \(m_i\) with \(H_\mathtt{eff,i} (\mathbf{m}, \mathbf{u})\). The Stratonovich type of the stochastic integral of the diffusion term ensures that each state process \(m_i\) takes values in \({{\mathbb {S}}}^2\). We refer to [3, 4] for further details on the model.

Our aim is to find a control process \(\mathbf{u}^*\) such that the solution process \(\mathbf{m}^*\) from (1.2) approximates a given deterministic profile \(\widetilde{\mathbf{m}} \equiv (\widetilde{m}_1, \ldots , \widetilde{m}_N) \in L^2\bigl ( 0,T; ({{\mathbb {S}}}^2)^N \bigr )\). More precisely, we aim to solve the following problem.

Problem 1.1

Let the parameters \(\delta , \nu , \alpha \ge 0\), and \( T, \lambda >0\), \(N \in {{\mathbb {N}}}\) as well as \(h \in C^2\bigl ( ({\mathbb S}^2)^N; {{\mathbb {R}}}\bigr )\) be given. Let \(\big (\Omega , \mathbb {P}, \mathcal {F},\{\mathcal {F}_t\}_{0 \le t \le T}\big )\) be a given stochastic basis with the usual conditions, and \(\mathbf{W}\) be a \(\{\mathcal {F}_t\}_{0 \le t \le T}\)-adapted \((\mathbb {R}^3)^N\)-valued Wiener process on it. Find a pair^{Footnote 1}

which minimizes the cost functional

subject to (1.2). We call such a minimizer a strong solution of the optimally controlled Landau–Lifschitz–Gilbert equation.

The problem (1.2) with \(\mathbf{u}=\mathbf{0}\) has been studied in [3]. In the case of a finite ensemble of nanomagnetic particles, the deterministic optimal control problem has been studied in [1, 2]; we also refer to [5] for a deterministic control problem for infinitely many ferromagnetic particles (a PDE). In [7], some of the present authors have studied Problem 1.1 in its weak form and constructed a weak optimal solution \(\Pi ^* := \bigl (\Omega ^*, {{\mathbb {P}}}^*, {\mathcal F}^*, \{{{\mathcal {F}}}^*_t\}_{0\le t\le T}, \mathbf{m}^*, \mathbf{u}^*, \mathbf{W}^*\bigr )\) of the underlying problem via the Young measure (relaxed control) approach for a compact control set \(\mathbb {U} \subset (\mathbb {R}^3)^N\) such that \(\mathbf{0}\in \mathbb {U}\), which may be generalized to the case \(\mathbb {U}=(\mathbb {R}^3)^N\) thanks to the coercivity of the cost functional with respect to \(\mathbf{u}\); see also [6] for an extension to infinite spin ensembles. To approximate it numerically, implementable strategies may be developed that rest on Pontryagin’s maximum principle which characterizes minimizers. In [7], a stochastic gradient method is proposed to generate a sequence of functional-decreasing approximate feedback controls, where the update requires to solve a coupled forward-backward SDE system. A relevant part here is to simulate a (time-discrete) backward SDE via the least-squares Monte-Carlo method, which requires significant data storage resources [6, 7], and thus limits the complexity of practically approachable Problem 1.1.

In this work, we use an alternative strategy which rests on the dynamic programming principle. This allows us to prove the existence of a unique strong solution of Problem 1.1, which sharpens results of [7]. Since the solution of the underlying SDE lies on \({{\mathcal {M}}} = ({{\mathbb {S}}}^2)^N\), the Hamilton–Jacobi–Bellman equation is defined on the manifold \([0,T] \times \mathcal {M}\). In Sect. 3 we verify that the value function is the unique solution of that Bellman equation and that it belongs to \(C^{1,2} \bigl ([0,T] \times \mathcal {M}\bigr )\). To solve this semi-linear PDE by deterministic numerical strategies seems non-accessible due to the high dimension of the underlying manifold \({{\mathcal {M}}} \); we also want to avoid a direct probabilistic representation of its solution which would involve a backward SDE. Indeed, we demonstrate how the nonlinear HJB equation may be replaced with a linear parabolic PDE (3.12) by applying the Hopf–Cole transformation. The quadratic form resp. linearity of the control in the cost functional resp. in the equation (1.2) together with the geometric character of the problem then lead to an isotropic quadratic term in the HJB equation (3.8), which is crucial for this transformation. The regularity of the value function and the optimal policy mapping is explicitly expressed through the regularity of the terminal condition h. Furthermore, the solution w of the linear parabolic PDE can now be represented via a Feynman–Kac formula. This is the starting point for the numerical scheme proposed in Sect. 5. To approximate the optimal pair \((\mathbf{m}^*, \mathbf{u}^*)\) numerically, a Monte-Carlo method for the solution w of the linear equation (3.8) and its tangential gradient \(\nabla _{\mathcal {M}}w\) is proposed, from which the optimal feedback function \(\bar{\mathbf{u}}\) can be obtained directly via (3.15). To approximate \(\nabla _{\mathcal {M}}w\) through a difference quotient with needed accuracy, we choose a stencil diameter \(\bar{h} = \mathcal {O}\big (1/\sqrt{M}\big )\) for a sufficiently large number of Monte-Carlo realizations M; see Remark 6.1. Importantly, this approach does not require larger data storage resources as [6, 7] does, but an ample calculation of iterates from related SDEs. Computational studies for the switching dynamics of single and multiple ferromagnetic particles are reported in Sect. 6.

The strategy which is proposed in this work to efficiently approximate the minimizer of the stochastic optimal control Problem 1.1 exploits its special structure to successfully apply the Hopf–Cole transformation—which then

1. (i)

   guarantees the existence and uniqueness of a classical solution of the associated Hamilton–Jacobi–Bellman equation, and of a strong solution of the optimal control problem,

2. (ii)

   leads to a characterisation of optimality through a linear parabolic PDE, which in turn allows the efficient numerical approximation in high-dimensional spaces through Monte-Carlo simulation without the need for computationally more costly backward SDEs.

The solution process \(\mathbf{m}\) of (1.2) attains values in \({{\mathcal {M}}} = (\mathbb {S}^2)^N\). For any

we have

For any \(\pmb {m} = ({m}_1, \ldots , {m}_N)^\top \in \mathcal {M}\), denote \(\pmb {\sigma }(\pmb {m}) = \mathtt{diag} \bigl (\sigma (m_1), \ldots , \sigma (m_N) \bigr )\), where \(\sigma ({m}_i) \in {\mathbb {R}}^{3 \times 3}\) is the matrix

Again, for any \((\pmb {m}, \pmb {u})\in \mathcal {M}\times (\mathbb {R}^3)^N\), we define

Then, using \(\pmb {\sigma }\) also here, and combining \(H_\mathtt{ani,i}({m}_i) + H_\mathtt{d,i}({m}_i)= - \mathtt{D}_i {m}_i\) with some \( \mathtt{D}_i= \mathtt{B}_i - \mathtt{A}\in {{\mathbb {R}}}^{3 \times 3}_\mathtt{diag} \) as well as \(\mathbf{D}=\mathtt{diag}\big ( \mathtt{D}_1, \ldots , \mathtt{D}_N\big )\), we have

where

The matrix \(\pmb {\Sigma }(\pmb {m}) \in {{\mathbb {R}}}^{3N \times 3N}\) is block-diagonal, with its i-th block

For \({m}_i \in \mathbb {S}^2\), one has

where \({{\mathcal {P}}}({m}_i)\) is the orthogonal projection onto the tangent plane of \(\mathbb {S}^2\) at \({m}_i\). Note that the diffusion term \(\nu \, m_i(s) \times (\circ \mathrm{d}W_i(s))\) in (1.2) can be re-written as \(\nu \, \sigma (m_i(s))\circ \mathrm{d}W_i(s)\). To state the dynamic programming equation, we introduce a family of stochastic Landau–Lifschitz–Gilbert equations with different initial times \(t \in [0,T]\) and states \(\pmb {m} \in {{\mathcal {M}}}\):

where \(\mathbf{f}\) is defined in (2.1). The solutions \(\mathbf{m} = \mathbf{m}^{t, {\pmb m}}\) of (2.3) thus depend on t and \({\pmb m}\); however, we shall drop the superscript of \(\mathbf{m}^{t, {\pmb m}}\) in the subsequent text for the ease of notation. For every \(0\le t\le T\), \(\mathbf{m}(t)=\pmb {m}\), and \(\mathbf{u} \in L^2_{\{ {{\mathcal {F}}}_s\}}\bigl ( \Omega ; L^2(t,T; (\mathbb {R}^3)^N)\bigr )\), there exists a unique strong solution \(\mathbf{m}=(m_1, \ldots , m_N) \in L^2_{\{\mathcal {F}_s\}}\big ( \Omega ; C(t,T; (\mathbb {R}^3)^N)\big )\) of (2.3). Indeed by considering the truncation of the control and then using the stochastic version of the Arzela-Ascoli theorem, Prohorov’s lemma and Jakubowski-Skorokhod representation theorem (cf. [14]), we have existence of a weak solution of (2.3). Moreover, an application of Gyöngy-Krylov’s characterization of convergence in probability introduced in [12] along with pathwise uniqueness of weak martingale solutions gives existence of a unique strong solution, see [6, Appendix]. Furthermore, by applying Itô’s formula to the functional \(\mathbf{x}\rightarrow \Vert \mathbf{x}\Vert _{\mathbb {R}^3}^2\) for \(m_i\) and using the vector identity \(\langle \mathbf{a}, \mathbf{a} \times \mathbf{b}\rangle =0\) for any \(\mathbf{a}, \mathbf{b}\in \mathbb {R}^3\), we have \(\mathbb {P}\)-a.s.,

Since \(m_i(t)\in \mathbb {S}^2\), we see that \(\mathbb {P}\)-a.s., each \(m_i\) is \(\mathbb {S}^2\)-valued, and thus \(\mathbf{m}\in L^2_{\{\mathcal {F}_s\}}\big ( \Omega ; C(t,T; (\mathbb {S}^2)^N)\big )\).

Because the paths of the Landau–Lifschitz–Gilbert process stay on the manifold \(\mathcal {M}\), the natural domain for the value function of the control problem is \([0,T] \times \mathcal {M}\). In order to make the connection between the controlled process \(\mathbf{m}\) on the one hand and the Hamilton–Jacobi–Bellman PDE posed on \([0,T] \times \mathcal {M}\) on the other hand, it is convenient to describe properties of \(\mathbf{m}\) purely in terms of quantities that are intrinsically defined on \(\mathcal {M}\), without referring to the ambient space \((\mathbb {R}^3)^N\). Of particular interest is Dynkin’s formula.

We begin by rewriting Itô’s formula with tangential derivatives. For \(1 \le \ell \le 3N\), let \(\pmb {\sigma }(\pmb {m})_\ell \) be the \(\ell \)th row of \(\pmb {\sigma }(\pmb { m})\), for \(\pmb {m} \in {\mathcal M}\). Then \(\partial _{\pmb {\sigma }(\pmb {m})_\ell }\) denotes the tangential derivative in the direction \(\pmb {\sigma }(\pmb {m})_\ell \), and \(\partial _{\pmb {\sigma }(\pmb {m})} := \bigl ( \partial _{\pmb {\sigma }(\pmb {m})_1}, \ldots , \partial _{\pmb {\sigma }(\pmb { m})_{3N}}\bigr ) \in [T_{\pmb {m}} {{\mathcal {M}}}]^{3N}\). Similarly, \(\partial _{\mathbf{f}(\pmb {m}, \pmb {u})}\) is the tangential derivative in the direction \(\mathbf{f}(\pmb {m}, \pmb {u})\).

We wish to apply Itô’s formula to \(\psi (s,\mathbf{m}(s))\) for any \(\psi \in C^{1,2}\big ([0,T]\times \mathcal {M}\big )\). One may directly return to the standard formula on \((\mathbb {R}^3)^N\), cf. [13, Chapter V.1], by extending \(\psi \) via

to \([0,T] \times (\mathbb {R}^3\setminus \{0\})^N\), for any \(\hat{\pmb {m}}=\big (\hat{{m}}_1,\ldots , \hat{{m}}_N\big )\in (\mathbb {R}^3\setminus \{0\})^N\). Then \({{\mathbb {P}}}\)-a.s.:

For \(\psi \in C^{1,2}([0,T]\times {{\mathcal {M}}})\), we associate the generator of the Markov process \(\mathbf{m}\) as

where \(\Delta _{{\mathcal {M}}}\) denotes the Laplace-Beltrami operator on \({{\mathcal {M}}}\), and define the operator

Lemma 2.1

For any \(\psi \in C^{1,2}\big ([0,T]\times \mathcal {M}\big )\) the process \(\psi (t, \pmb {m})\) satisfies Dynkin’s formula

Proof

Re-writing (2.4) with \(s=T\) in Itô form, we have

As in the proof of [18, Proposition 3.2] we conclude that

Taking the expectation then leads to (2.6), recalling that the Itô integral is a martingale. \(\square \)

For any \((t,\pmb {m})\in [0,T]\times \mathcal {M}\), we consider problem (2.3) to now construct the associated Hamilton–Jacobi–Bellman equation, following the formal rules of dynamic programming. We then use the Hopf–Cole transformation to replace the nonlinear HJB equation by a linear PDE and show the existence of a unique classical solution, which then implies existence of a unique classical solution of the original nonlinear HJB equation. Next, we present a verification theorem which shows that the value function is indeed equal to the solution of the HJB equation. We describe the optimal control through an optimal feedback function which is written explicitly in terms of the value function.

Let us define the Lagrangian

where the parameters \(\delta , \lambda \) are given in Problem 1.1, and \(L\big (\pmb {m}, \pmb {u}\big )\) appears in the cost functional.

Let \(\big (\Omega , \mathbb {P}, \mathcal {F},\{\mathcal {F}_s\}_{t \le s \le T}\big )\) be a given filtered probability space satisfying the usual hypotheses, and \(\mathbf{W}\) is a \(\{\mathcal {F}_s\}_{t \le s \le T}\)-adapted \((\mathbb {R}^3)^N\)-valued Wiener process on it. We denote by \(\mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) the set of all admissible pairs \((\mathbf{m}, \mathbf{u})\) such that \(\mathbf{u} \in L^2_{\{{{\mathcal {F}}}_{s}\}}\bigl ( \Omega ; L^2( t,T; ({\mathbb R}^{3})^N )\bigr )\), and \(\mathbf{m}{(\cdot )}\) is the unique \(\{ \mathcal {F}_s \}_{t \le s \le T}\)-adapted \( \mathcal {M}\)-valued strong solution of (2.3). In fact, the superscript \(\pmb {s}\) refers to the fact that we search an optimal admissible pair on the given filtered probability space. It follows for admissible \((\mathbf{m}, \mathbf{u})\) that \(L\big (\mathbf{m}(\cdot ), \mathbf{u}(\cdot )\big )\in L^1_{\{{{\mathcal {F}}}_{s}\}}\bigl ( \Omega ; L^1( t,T; \mathbb {R})\bigr )\) and \(h(\mathbf{m}(T))\in L^1_{\mathcal {F}_T}(\Omega ; \mathbb {R})\), recalling that h is a continuous function on a compact manifold. As a further consequence, Dynkin’s formula (2.6) then holds for all \(\psi \in C^{1,2}([t,T]\times \mathcal {M})\) and \((\mathbf{m}, \mathbf{u}) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\).

Our aim is to achieve \(\mathbf{m}\) close to a reference configuration \(\widetilde{\mathbf{m}} \in C^2([t,T]; \mathcal {M})\) by selecting an optimal control \(\mathbf{u}^*\). The cost functional on \(\mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) is

We write the value function of \(\mathcal {J}\) as

Note that, thanks to [3, Proposition 1.33], there exists a unique strong solution \(\mathbf{m}\) of (2.3) for \(\mathbf{u}=\mathbf{0}\), and hence the value function is uniformly bounded since h is a given continuous function and \(\Vert \mathbf{m}\Vert _{(\mathbb {R}^3)^N}^2=N\).

3.1 The Hamilton–Jacobi–Bellman Equation

We define the Hamiltonian

Note that in the definition of \(\mathcal {H}\), the letter \(\pmb {m}\) stands for the identity map \(\pmb {m} \mapsto \pmb {m}\). Using (2.1), we evaluate the supremum analytically. We have for the tangential gradient \(\nabla _\mathcal {M} \psi \)

where \(^*\) denotes the convex conjugate function, and \(\mathbf{Q}\) is the \(3N \times 3N\) matrix given by

Since \(\nabla _\mathcal {M} \psi (\pmb {m})\) belongs to the tangent space \(T_{\pmb {m}}\mathcal {M}\), it follows that

Note that for any \(\pmb {m} \in \mathcal {M}\), \(\Vert \pmb {m} \times \nabla _\mathcal {M} \psi (\pmb {m}) \Vert _{(\mathbb {R}^3)^N} = \Vert \nabla _\mathcal {M} \psi (\pmb {m}) \Vert _{(\mathbb {R}^3)^N}\). Thus by Pythagoras’ theorem, we have

Let us denote

Then, in summary, we have

We point out that (3.6) ensures that the quadratic term in the Hamiltonian is isotropic, which is crucial for the Hopf–Cole transformation below. We now state the Hamilton–Jacobi–Bellman equation, whose solution we denote by W:

The Hopf–Cole transformation: The HJB equation (3.8) is a second-order nonlinear PDE on a high-dimensional domain and therefore without further understanding of its structure computationally expensive to solve; the study of its wellposedness as well as the regularity of its solution is non-trivial. We use the Hopf–Cole transformation \(w = \exp \bigl (- \beta \, W \bigr )\), \(\beta \in \mathbb {R}\) given below, to substitute (3.8) by the linear PDE (3.12).

We span the tangent space of \(\mathcal {M}\) at any point \(\pmb {m}\) by the orthonormal tangent vectors

It is convenient to conceptually let \(\partial _{i,1}, \partial _{i,2}\) span the local coordinates associated to the ith sphere. Then we have the following relations: for \(j\in \{1,2\}\) and \(i\in \{1,2,\cdots , N\}\),

We see that

Therefore, we have the following correspondences:

Recalling \(\Vert \nabla _\mathcal {M} W \Vert _{(\mathbb {R}^3)^N}^2 = \sum _i \bigl [ |\partial _{i,1} W|^2 + |\partial _{i,2} W|^2 \bigr ]\), we choose \( \displaystyle \beta = \frac{ C_\mathtt{ext}^2(1 + \alpha ^2)}{\lambda \nu ^2}\) to obtain a cancellation of the quadratic nonlinearity via (3.10). Substituting the respective terms in (3.8) and multiplying by \( - \beta \, w \ne 0\), we arrive at the second-order linear equation

where \(c(t,\pmb {m})=\beta \,\delta \, \Vert \pmb {m} - \widetilde{\mathbf{m}}(t) \Vert _{(\mathbb {R}^3)^N}^2\). The following theorem shows that weak solutions w may be examined in the Sobolev-Bochner space

A weak solution of (3.12) is a \(w \in W^{1,2,2}\big ([0,T]; H^1({\mathcal {M}}), H^{-1}({\mathcal {M}})\big )\) such that for all \(\psi \in L^2\big (0,T; H^1(\mathcal {M})\big )\):

Theorem 3.1

There exists a unique classical solution \(w \in C^{1,2}\big ([0,T] \times \mathcal {M}\big )\) of (3.12) which is also the unique weak solution in \(W^{1,2,2}\big ([0,T]; H^1({\mathcal {M}}), H^{-1}({\mathcal {M}})\big )\). The function

is the unique classical solution of the Bellman equation (3.8).

Proof

The existence of a unique solution w of (3.13) in \(W^{1,2,2}\big ([0,T]; H^1({\mathcal {M}}), H^{-1}({\mathcal {M}})\big )\) is for instance given in [19, p. 224]. Using charts of \(\mathcal {M}\) and a partition of unity to represent (3.13) locally on the flat space, and an application of parabolic regularity theory ensures that w belongs to \(C^{1,2}\big ([0,T] \times \mathcal {M}\big )\).

We shall now establish that w is positive so that (3.14) may be applied. Let \(\phi (t,x) = \exp (\gamma t) [w(t,x) - \epsilon ]\) with \(\epsilon := \min _{\tilde{x}} w(T,\tilde{x}) > 0\) and \(\gamma := \Vert \mathrm{div} \mathbf{b} \Vert _\infty + 1\). In (3.13) we substitute \(\psi \) by \(\exp (\gamma t) \psi \) and then add and subtract \(\gamma \phi \cdot \psi \). Applying the product rule one obtains

Now let \(\psi = \min \{0, \phi \}\). Where \(\psi \ne 0\) one has \(\phi = \psi \) so that

Thus

By construction \(\psi (T) \equiv 0\). We conclude that \(\psi \equiv 0\) and therefore \(\phi \ge 0\). Thus \(w \ge \epsilon > 0\).

The above implies that the nonlinear HJB equation (3.8) has a classical solution W. Moreover, the construction of the Hopf–Cole transformation directly ensures that a function \(w \in C^{1,2}\big ([0,T] \times \mathcal {M}\big )\) is a classical solution of (3.12) if and only if W of (3.14) solves (3.8) classically, which then implies the existence of a unique classical solution W, given by (3.14), of the HJB equation (3.8). \(\square \)

It is easy to see that additional smoothness of the terminal condition h directly translates into additional regularity of w and W.

3.2 The Verification Theorem

In this subsection, we show that \(V = W\), i.e., the value function of the optimal control problem is equal to the solution of the above HJB equation (3.8). The following verification theorem also provides an explicit formula for the optimal control, which inherits the smoothness of \(\nabla _{\mathcal {M}} W\).

Theorem 3.2

The value function \(V(t,\pmb {m})\) in (3.3) is the unique classical solution of the nonlinear HJB equation (3.8):

Problem 1.1 admits a minimizer \((\mathbf{m}^*, \mathbf{u}^*)\in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) such that \(\mathcal {J}\big (t, \pmb {m};(\mathbf{m}^*, \mathbf{u}^*)\big )=V(t,\pmb {m})\) and \(\mathbf{u}^*(s)=\bar{\mathbf{u}}\big (s,\mathbf{m}^*(s)\big )\), where

Proof

The proof is divided into two steps.

Step 1: First we show that \(W \le V\). Let \((\mathbf{m}, \mathbf{u})\in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) be any admissible pair. Now, for any \(\psi \in C^{1,2}\big ([0,T) \times \mathcal {M}\big )\), we have, thanks to the definition of the Hamiltonian \(\mathcal {H}\), (3.1), and (2.5),

and therefore one has

Because W is a smooth classical solution, we may substitute \(\psi \) by W, in which case the left-hand side of (3.17) vanishes. In other words,

Indeed, the existence of a classical solution W avoids a more complicated construction to arrive at a bound like (3.18), which would be necessary in a setting with viscosity solutions. Applying now Dynkin’s formula (2.6) with W in place of \(\psi \) and recalling the final time conditions, we conclude from (3.18),

Because this result holds for all admissible pairs \( (\mathbf{m},\mathbf{u}) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\), it follows that \(W \le V\).

Step 2: Now we show that there exists an admissible pair \( (\mathbf{m}^*, \mathbf{u}^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) such that \(W(t,\pmb {m})= \mathcal {J}\big (t,\,\pmb {m};(\mathbf{m}^*, \mathbf{u}^*) \big )\). Recalling (3.4), we find that the supremum in the definition of the Hamilton–Jacobi–Bellman operator is attained by the control

Since W is a \(C^{1,2}\)-solution of (3.8), we see that \(\bar{\mathbf{u}}\) is continuously differentiable and bounded on \([0,T] \times \mathcal {M}\). Moreover, \(\mathcal {M}\ni {\pmb m}\mapsto \bar{\mathbf{u}}(t,\pmb {m})\) is Lipschitz. Indeed, for any \(\pmb {m}_1,\, \pmb {m}_2 \in \mathcal {M}\)

Since \(W\in C^{1,2}\big ( [0,T]\times \mathcal {M}\big )\), by the mean-value theorem, either applied in combination with Whitney’s extension theorem [8, Section 6.5] in the ambient space \((\mathbb {R}^3)^N\) or directly to the geodesics of \(\mathcal {M}\), we have

Thus,

Now, on the given stochastic basis \((\Omega , \mathbb {P}, \mathcal {F},\{\mathcal {F}_s\}_{t \le s \le T})\) and for the \(\{\mathcal {F}_s\}_{t \le s \le T}\)-adapted Brownian motion \(\mathbf{W}(s)\), the stochastic differential equation

has a pathwise unique \(\mathcal {M}\)-valued solution, recall the argument for (2.3). Then the process

belongs to \( L^2_{\{{{\mathcal {F}}}_{s}\}}\bigl ( \Omega ; L^2( t,T; (\mathbb {R}^3)^N)\bigr )\). Thus \((\mathbf{m}^*, \mathbf{u}^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\). With this admissible pair \((\mathbf{m}^*, \mathbf{u}^*)\), the inequality in (3.16) turns into equality. Again, by using Dynkin’s formula along with initial data \(\mathbf{m}^*(t)=\pmb {m}\), we see that

Recalling that by Theorem 3.1 the HJB equation (3.8) has a unique solution W, we conclude from the above that \(V = W\). \(\square \)

Now we show the uniqueness of the optimal admissible pair \((\mathbf{m}^*, \mathbf{u}^*)\), and thus in particular of the strong solution of Problem 1.1. We remark that the uniqueness of the optimal admissible pair is not automatically provided from the uniqueness of solutions of the HJB equation, which instead was used to characterize the value function through the differential equation (2.3).

Theorem 3.3

The pair \((\mathbf{m}^*, \mathbf{u}^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) constructed in Theorem 3.2 is the unique minimizer of \(\mathcal {J}\big (t,\pmb {m}; (\cdot ,\cdot )\big )\) in the sense that if there exists any other optimal pair \((\mathbf{m}_1^*, \mathbf{u}_1^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\), then \( \mathbf{m}^*(s)= \mathbf{m}_1^*(s)\) and \(\mathbf{u}^*(s)=\mathbf{u}_1^*(s)\) for a.e. \(s\in [t,T],\,\mathbb {P}\)-a.s.

Proof

Step 1: Let \((\mathbf{m}^*, \mathbf{u}^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) be an optimal pair. Similar to [20, Chapter 5, Theorem 5.1], we note that then (3.19) holds as equality with \((\mathbf{m}^*, \mathbf{u}^*)\) in place of \((\mathbf{m}, \mathbf{u})\). This implies that also (3.16) holds with equality for a.e. \(s\in [t,T]\), \(\mathbb {P}\)-a.s. Hence \((\mathbf{m}^*, \mathbf{u}^*)\) and thus every optimal pair satisfies

for a.e. \(s\in [t,T]\), \(\mathbb {P}\)-a.s.

Step 2: Suppose there exists another optimal pair \((\mathbf{m}_1^*, \mathbf{u}_1^*) \in \mathcal {U}_{\pmb {m}}^{\pmb {s}}[t,T]\) and let \(\widetilde{\mathbf{m}}^*= \mathbf{m}^* - \mathbf{m}_1^*\). Then \(\widetilde{\mathbf{m}}^*\) is a strong solution to the following SDE: for \(s\in (t,T]\)

with \(\widetilde{\mathbf{m}}^*(t)=\mathbf{0}\). In view of (2.2), we observe that

Thus, the Eq. (3.21) reduces to

We now apply Itô’s formula to the functional \(\mathbf{x}\mapsto \Vert \mathbf{x}\Vert _{(\mathbb {R}^3)^N}^2\) for the above equation, and then use (3.20) to have

for some constant \(C>0\). An application of Gronwall’s lemma then implies \(\mathbb {P}\)-a.s., \(\mathbf{m}^*(s)=\mathbf{m}_1^*(s)\), and therefore from (3.20) we get \(\mathbf{u}^*(s)=\mathbf{u}_1^*(s)\) for a.e. \(s\in [t,T]\). Thus, the optimal control problem admits a unique strong solution, which is an improvement over [7]. \(\square \)

Remark 3.4

In [7], \(\mathbb {P}\)-a.s. the orthogonality of an optimal state and control was shown both theoretically and numerically. In our approach, we also see \(\mathbb {P}\)-a.s. the orthogonality of \(\mathbf{m}^*\) and \(\mathbf{u}^*\). Indeed, by using the vector identity \(\langle \mathbf{a}, \mathbf{a} \times \mathbf{b}\rangle =0\) for any \(\mathbf{a}, \mathbf{b}\in (\mathbb {R}^3)^N\), and the fact that \({\pmb m}\in \mathcal {M}\) and \(\nabla _{\mathcal {M}}W(\cdot ,\pmb {m})\) are orthogonal, we have \(\mathbb {P}\)-a.s., from (3.15)

For its computational evidence, see Fig. 3 for an ensemble of \(N=3\) particles.

To solve the linear PDE (3.12) numerically by deterministic methods is still demanding since it is posed on \(\mathcal {M}\subset (\mathbb {R}^3)^N\). Therefore, we choose a probabilistic representation of the solution of (3.12) which requires to solve the following forward stochastic differential equation, defined on a given stochastic basis \((\Omega , \mathbb {P}, \mathcal {F},\{\mathcal {F}_s\}_{t \le s \le T})\) with 3N-dimensional Brownian motion \(\mathbf{W}\),

where \(\mathbf{b}(\cdot )\) is defined in (3.7). Equation (4.1) has a strong solution \(\pmb {{\mathfrak {m}}}\) taking values in \(\mathcal {M}\). Let \(G(s)= \displaystyle \exp \Big (-\int _t^s R\big (r, \pmb {{\mathfrak {m}}}(r)\big )\,\mathrm{d}r\Big )\) where \(\displaystyle R(r,\pmb {{\mathfrak {m}}}(r))= \beta \,\delta \Vert \pmb {{\mathfrak {m}}}(r) - \widetilde{\mathbf{m}}(r) \Vert _{(\mathbb {R}^3)^N}^2\). By using the Itô product rule applied to \(G(s)\,w\big (s, \pmb {{\mathfrak {m}}}(s)\big )\), where w is the classical solution of the linear parabolic PDE (3.12), we arrive at the following Feynman–Kac representation [15, Theorem 7.6] for the solution of (3.12) with terminal datum \( \exp ( -\beta \, h(\pmb {m})) = \exp ( -C_\mathtt{ext}^2(1 + \alpha ^2) (\lambda \nu ^2)^{-1} \, h(\pmb {m}))\):

We note that because of the linearity of (3.12) the Feynman–Kac representation can be used in place of a backward SDE.

4.1 A Numerical Scheme for (4.1)

To approximate the solution \(\pmb {{\mathfrak {m}}}\) of (4.1), we use the semi-implicit method proposed in [17]. Now \(\mathbf{b}(\pmb {m})\) from (3.7) can be re-written as

where \(\mathbf{Q}\) is defined in (3.5). Let \(T>0\) be fixed. For \(J\in \mathbb {N}\), let \(I_J^0:= \{ t_j\}_{j=0}^J\) be a partition of [0, T] with time step size \( \displaystyle \tau = {T/J} >0\). Let \(I_J^{\ell } \subset I_J^0\) be the sub-partition on \([t_{\ell },T]\), where \(\ell \in \{0,1,\ldots , J-1 \}\). Let \(\{\pmb {\xi }^j\}_{j=\ell }^J\) is a \(({{\mathbb {R}}}^{3})^N\)-valued random walk of \(I_J^{\ell }\) with \(\pmb {\xi }^j := (\xi ^j_1, \ldots , \xi ^j_{N})\), where (\(1\le i\le N\)) \(\xi ^j_i=\big ( \xi ^j_{i,l}\big )_{1\le l\le 3}\) are i.i.d. \({\mathbb R}^3\)-valued (discrete) random variables such that each

1. i)

   \({\xi }^j_{i,l}\) satisfies \(\mathbb {E}\big [{\xi }^j_{i,l}\big ]=0\) and \(\mathbb {E}\big [\big |{\xi }^j_{i,l}\big |^2\big |\big ]=\tau \),

2. ii)

   for every integer \(p\ge 1\), there exists \(C_p>0\) such that \(\mathbb {E}\big [ |{\xi }^j_{i,l}|^{2p}\big ] \le C_p \tau ^p\).

Let now \(\ell \) be fixed and \(\pmb {{\mathfrak {m}}}^{\ell }= \pmb {m}\in \mathcal {M}\) be given. We determine the \(\mathcal {M}\)-valued random variables \(\{\pmb {{\mathfrak {m}}}^j \}_{j=\ell }^J\) via \(\big (i=1,2,\ldots , N\big )\)

where \(\pmb {{\mathfrak {e}}}^j=\big ({{\mathfrak {e}}}_1^j, {\mathfrak e}_2^j,\ldots , {{\mathfrak {e}}}_N^j\big )\), and the function \(\bar{\mathtt{a}}_i\) for each \(i=1,2,\ldots , N\) is given by

Note that (4.3a)–(4.3b) is a system of linear equations, leading to fast simulation times. Furthermore, the numerical schemes ensures \(\pmb {{\mathfrak {m}}}^j\) takes values on \(\mathcal {M}\). This is exploited when applying arguments from [3] in this finite ensemble setting to conclude weak convergence of iterates \(\{\pmb {{\mathfrak {m}}}^j \}_{j=\ell +1}^J\), which here are driven by a random walk, towards the strong solution of (4.1) for \(\tau \rightarrow 0\).

In order to simulate the optimal pair \((\mathbf{m^*}, \mathbf{u}^*)\), we need to solve the equations (4.1), (3.15) and (2.3) numerically. To solve (2.3) numerically, we replace \(\mathbf{u}\) by \(\bar{\mathbf{u}}(s,\cdot )\) defined in (3.15).

5.1 HJB Solution

The classical solution w of (3.12) is given by (4.2). In order to approximate it, and hence the classical solution W of the nonlinear HJB equation (3.8), we proceed as follows:

1. a)

   Compute all the iterates \(\{\pmb {{\mathfrak {m}}}^j\}_{j=\ell }^J\) via (4.3a)–(4.3b) along \(I_J^{\ell }\) and store them.

2. b)

   Approximate the integral in (4.2) by Gauss-Legendre quadrature [16, Section 10.3], where we use the piecewise affine interpolation of the iterates \(\{\pmb {{\mathfrak {m}}}^j\}_{j=\ell }^J\) via

   $$\begin{aligned} \pmb {{\mathfrak {m}}}(r):= \frac{r-t_j}{\tau } \pmb {{\mathfrak {m}}}^{j+1} + \frac{t_{j+1}-r}{\tau }\pmb {{\mathfrak {m}}}^j \qquad \big (r\in [t_j, t_{j+1})\big )\,. \end{aligned}$$

   (5.1)
3. c)

   Since \(\widetilde{\mathbf{m}}\in C^2([t_{\ell },T]; \mathcal {M})\), we use the piecewise affine interpolation \(\widetilde{\mathbf{m}}(r)\) of the iterates \(\widetilde{\mathbf{m}}^j\equiv \widetilde{\mathbf{m}}(t_j)\).

4. d)

   Approximate \(\mathbb {E}_{t,\pmb {m}}\) in (4.2) via Monte-Carlo estimation along with the variance reduction method of antithetic variates (see e.g. [11, Subsection 4.2]).

Thus, we can simulate the quantities \(w(t_{\ell }, \pmb {\mathfrak m}^\ell )\), and hence \(W(t_{\ell }, \pmb {{\mathfrak {m}}}^\ell )= \displaystyle -\frac{1}{\beta } \log (w(t_{\ell }, \pmb {\mathfrak m}^\ell ))\).

5.2 Optimal Feedback Transformation

To approximate the function \(\bar{\mathbf{u}}\) at any point \((t,\pmb {m})\), we need to approximate \(\nabla _\mathcal {M} W(t,\pmb {m})\), which again demands to approximate \(\nabla _\mathcal {M} w(t,\pmb {m})\) thanks to the Hopf–Cole transformation. For the latter, we may proceed in two different ways:

1. i)

   Method A: We take the expectation first and then use the central difference quotient to approximate the tangential gradient. More precisely, for any \(\pmb {m}=\big ({m}_1,{m}_2,\ldots , {m}_N\big )\in \mathcal {M}\) and \(\bar{h}>0\), define for \(i=1,2,\ldots , N\),

   $$\begin{aligned} \begin{aligned} \pmb {m}_{\bar{h},i,l}^+: = \Big ( {m}_1, \ldots , \frac{{m}_i + \bar{h} \mathbf{e}_l}{\Vert {m}_i + \bar{h} \mathbf{e}_l\Vert _{\mathbb {R}^3}}, \ldots , {m}_N\Big )\,, \\ \pmb {m}_{\bar{h},i,l}^-: = \Big ( {m}_1, \ldots , \frac{{m}_i - \bar{h} \mathbf{e}_l}{\Vert {m}_i - \bar{h}{} \mathbf{e}_l\Vert _{\mathbb {R}^3}}, \ldots , {m}_N\Big )\,, \end{aligned} \end{aligned}$$

   (5.2)

   where \(\mathbf{e}_l\) is the l-th identity vector in \(\mathbb {R}^3,\, 1\le l\le 3\). Recall that, by using Sect. 5.1, we can calculate \(w(t, \pmb {m}_{\bar{h},i,l}^+)\) and \(w(t, \pmb {m}_{\bar{h},i,l}^-)\) for any \(t\in [0,T]\); therefore, we approximate \(\nabla _\mathcal {M} w(t,\pmb {m})\) by \( \nabla _\mathcal {M} w(t,\pmb {m})\simeq \Big ( D_1 w(t,\pmb {m}),\, \ldots ,\, D_N w(t,\pmb {m})\Big )\), where for any \(i=1,2,\ldots , N\), and \(l=1,2,3\),

   $$\begin{aligned}&D_i w(t,\pmb {m}):= \Big (d_{1,i} w(t, \pmb {m}),\, d_{2,i} w(t, \pmb {m}),\, d_{3,i}w(t, \pmb {m})\Big )^\top \,, \\&d_{l,i} w(t,\pmb {m}):= \frac{1}{2\bar{h}}\Big [ w(t,\pmb {m}_{\bar{h},i,l}^+ )- w(t, \pmb {m}_{\bar{h},i,l}^-)\Big ]\,. \end{aligned}$$

   Hence, we approximate \(\nabla _{\mathcal {M}}W(t_{\ell },\pmb {{\mathfrak {m}}}^\ell )\) by

   $$\begin{aligned} \nabla _\mathcal {M} W(t_{\ell },\pmb {{\mathfrak {m}}}^{\ell })\simeq -\frac{1}{\beta w(t_{\ell },\pmb {{\mathfrak {m}}}^{\ell })} \Big ( D_1 w(t_{\ell },\pmb {{\mathfrak {m}}}^{\ell }),\, \ldots ,\, D_N w(t_{\ell },\pmb {{\mathfrak {m}}}^{\ell })\Big )\,. \end{aligned}$$

   (5.3)
2. ii)

   Method B: In contrast to Method A, we first use the central difference quotient and then take the expectation to approximate the gradient \(\nabla _\mathcal {M} w\). For all \((t, \pmb {m})\in [0,T)\times \mathcal {M}\), we define the random variable

   $$\begin{aligned} H(t,{\pmb {m}}):= \exp \big ( -\beta \, h(\pmb {{\mathfrak {m}}}(T))\big )\exp \Big (- \beta \,\delta \int _t^T \Vert \pmb {{\mathfrak {m}}}(r)-\widetilde{\mathbf{m}}(r) \Vert _{(\mathbb {R}^3)^N}^2\,\mathrm{d}r\Big )\,, \end{aligned}$$

   (5.4)

   where \(\pmb {{\mathfrak {m}}}\) solves (4.1) with \(\pmb {{\mathfrak {m}}}(t)=\pmb {m}\). Let \(\pmb {m}_{\bar{h},i,l}^+\) and \(\pmb {m}_{\bar{h},i,l}^-\) be the points in \(\mathcal {M}\) as defined above. We compute the central difference quotients component-wise, and use (4.3a)–(4.3b) to approximate related solutions from (4.1) in (5.4). For \(i=1,2,\ldots , N\), and \(l=1,2,3\), we define

   $$\begin{aligned} \begin{aligned} d_{l,i} H(t,\pmb {m})&:= \frac{1}{2\bar{h}}\Big [ H(t,\pmb {m}_{\bar{h},i,l}^+ )- H(t, \pmb {m}_{\bar{h},i,l}^-)\Big ] \,, \\ d_i H(t,\pmb {m})&:= \Big (\mathbb {E}\big [ d_{1,i}H(t, \pmb {m})\big ],\, \mathbb {E}\big [ d_{2,i}H(t, \pmb {m})\big ],\, \mathbb {E}\big [ d_{3,i}H(t, \pmb {m})\big ]\Big )^\top \,. \end{aligned} \end{aligned}$$

   (5.5)

   We approximate the expectation in (5.5) via Monte-Carlo estimation together with the method of antithetic variates. We then approximate \(\nabla _\mathcal {M} w(t,\pmb {m})\) (hence \(\nabla _\mathcal {M} W(t,\pmb {m})\)) as

   $$\begin{aligned}&\nabla _\mathcal {M} w(t,\pmb {m})\simeq \Big ( d_1 H(t,\pmb {m}),\, \ldots ,\, d_N H(t,\pmb {m})\Big )^\top \,, \\&\nabla _\mathcal {M} W(t,\pmb {m})\simeq -\frac{1}{\beta w(t,\pmb {m})}\Big ( d_1 H(t,\pmb {m}),\, \ldots ,\, d_N H(t,\pmb {m})\Big )^\top \,. \end{aligned}$$

Thus, we simulate the transformation function \(\bar{\mathbf{u}}(t_{\ell }, \pmb {{\mathfrak {m}}}^\ell )\) from (3.15) by using one of the above methods, where the sequence \((t_{\ell }, \pmb {{\mathfrak {m}}}^\ell )\) is described in Sect. 4.1. Therefore, while the approximation of w depends only on the numerical parameter M, for \(\nabla _\mathcal {M} w\) the approximation depends on the numerical parameters M and \(\bar{h}\).

5.3 Optimal State

We use again the semi-implicit method proposed in [17] to approximate the solution in (2.3) in which each realization takes values in \({{\mathcal {M}}}\). For any \(\pmb {m}=({m}_1,\ldots , {m}_N)\in \mathcal {M}\), and \(i=1,2,\ldots ,N\), define

where \(\mathbf{Q}\) is defined in (3.5). We use again the scheme (4.3a)–(4.3b), and Sect. 5.2 to find a \(\mathcal {M}\)-valued random variables \(\{\mathbf{m}^j \}_{j=\ell }^J\) along \(I_J^{\ell }\) with \(\mathbf{m}^{\ell }= \pmb {m}\in \mathcal {M}\) and \(\ell \in \{0,1,\ldots , J-1 \}\), where \(\bar{\mathtt{a}}_i\big (\pmb {{\mathfrak {m}}}^j\big )\) resp. \(\bar{\mathtt{a}}_i\big (\frac{\pmb {{\mathfrak {m}}}^j + \pmb {{\mathfrak {e}}}^j}{2}\big )\) in (4.3a) resp. (4.3b) is replaced by \(\mathrm{a}_i\big (\mathbf{m}^j, \bar{\mathbf{u}}(t_j, \mathbf{m}^j)\big )\) resp.

such that the iterates \(\{\mathbf{m}^j\}_{j=\ell }^J\) converge towards a weak solution of (2.3) for \(\tau \rightarrow 0\). Moreover, the iterates \(\{\mathbf{u}^j=\bar{\mathbf{u}}(t_j, \mathbf{m}^j)\}_{j={\ell }}^J\) defines the discrete optimal control along \(I_J^\ell \).

In summary, we have the following algorithm to compute the optimal solution and control along with Method B.

Algorithm 5.1

Let \(\mathbf{m}^{0} \in (\mathbb {S}^{2})^{N}\), \(T > 0\), \(M \in \mathbb {N}\) be given. For \(J \in \mathbb {N}\), let \(I_{J}^{0} := \{ t_{j} \}_{j=0}^{J}\) be a partition of [0, T] with time step size \(\tau = \frac{T}{J} > 0\). Denote by \(I_{J}^{\ell } \subset I_{J}^{0}\) the sub-partition on \([t_{\ell },T]\), where \(\ell \in \{ 0,1,\ldots ,J-1 \}\). Let now \(\ell \) be fixed and \(\mathbf{m}^{\ell } = \pmb {m} \in \mathcal {M}\) be given.

1. (I)

   Compute M-samples \(\mathcal {S}_{\pmb {\xi }}^{\ell } := \{ \mathcal {S}_{\pmb {\xi }}^{\ell ,k} \}_{k=1}^{M}\), \(\mathcal {S}_{\pmb {\xi }}^{\ell ,k} := \{ \pmb {\xi }^{j}(\omega _{k}) \}_{j=\ell }^{J}\) on \(I_{J}^{\ell }\).

2. (II)

   For \(i = 1,\ldots ,N\) do:

   For \(l = 1,2,3\) do:

   1. (1)

      Based on \(\pmb {m}:=\mathbf{m}^{\ell }\), compute \(\pmb {m}_{\bar{h},i,l}^{+}\) and \(\pmb {m}_{\bar{h},i,l}^{-}\) as in (5.2).

   2. (2)

      For \(k = 1,\ldots ,M\) do:

      1. (a)

         Compute \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{+,\ell ,k} := \{\pmb {{\mathfrak {m}}}^{j}(\omega _{k})\}_{j=\ell }^{J}\) resp.  \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{-,\ell ,k} := \{\pmb { \mathfrak m}^{j}(\omega _{k}) \}_{j=\ell }^{J}\) on \(I_{J}^{\ell }\) via scheme (4.3a)–(4.3b) for \(\bar{\texttt {a}}_{i}(\pmb {{\mathfrak {m}}}^{j})\) and \(\bar{\mathtt{a}}_i\big (\frac{\pmb {{\mathfrak {m}}}^j + \pmb {{\mathfrak {e}}}^j}{2}\big )\) using \(\{ +\pmb {\xi }^{j}(\omega _{k}) \}_{j=\ell }^{J}\) resp.  \(\{ -\pmb {\xi }^{j}(\omega _{k}) \}_{j=\ell }^{J}\).

      2. (b)

         Compute \(H(t_{\ell }, \pmb {m}_{\bar{h},i,l}^{+},\omega _{k})\) resp.  \(H(t_{\ell },\pmb {m}_{\bar{h},i,l}^{-},\omega _{k})\) in (5.4) based on \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{+,\ell ,k}\) resp.  \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{-,\ell ,k}\) to determine \(d_{l,i}H(t_{\ell },\pmb {m},\omega _{k})\) in (5.5) and store it.

   3. (3)

      Approximate \(\mathbb {E}[d_{l,i}H(t_{\ell },\pmb {m})]\) in (5.5) via Monte-Carlo estimation along with the variance reduction method of antithetic variates.

3. (III)

   Set \(\nabla _{\mathcal {M}} w(t_{\ell },\pmb {m}) \approx \bigl ( d_{1}H(t_{\ell },\pmb {m}),\ldots ,d_{N}H(t_{\ell },\pmb {m}) \bigr )^{\top }\) with \(d_{i}H(t_{\ell },\pmb {m})\) as in (5.5), and compute \(\bar{\mathbf{u}}(t_{\ell },\mathbf{m}^\ell )\) as in (3.15).

4. (IV)

   Compute \(\mathbf{e}^\ell \in (\mathbb {R}^{3})^{N}\) via Scheme (4.3a) with \(\mathrm{a}_{i}\bigl ( \mathbf{m}^{\ell }, \bar{\mathbf{u}}(t_{\ell },\mathbf{m}^{\ell }) \bigr )\) defined in (5.6).

5. (V)

   Define \(\mathbf{g}^{\ell } := \left( \frac{g_{1}^{\ell }}{{\Vert g_{1}^{\ell } \Vert }_{\mathbb {R}^{3}}},\ldots , \frac{g_{N}^{\ell }}{{\Vert g_{N}^{\ell } \Vert }_{\mathbb {R}^{3}}} \right) ^{\top }\) with \(g_{i}^{\ell } := \frac{m_{i}^{\ell }+e_{i}^{\ell }}{2}\,\,(1\le i\le N)\).

6. (VI)

   For \(i = 1,\ldots ,N\) do:

   For \(l= 1,2,3\) do:

   1. (4)

      Based on \(\pmb {g}:=\mathbf{g}^{\ell }\), compute \(\pmb {g}_{\bar{h},i,l}^{+}\) and \(\pmb {g}_{\bar{h},i,l}^{-}\) as in (5.2).

   2. (5)

      For \(k = 1,\ldots ,M\) do:

      1. (c)

         Compute \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{+,\ell ,k} := \{\pmb {{\mathfrak {m}}}^{j}(\omega _{k})\}_{j=\ell }^{J}\) resp.  \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{-,\ell ,k}:= \{\pmb { {\mathfrak {m}}}^{j}(\omega _{k})\}_{j=\ell }^{J}\) on \(I_{J}^{\ell }\) via scheme (4.3a)–(4.3b) for \(\bar{\texttt {a}}_{i}(\pmb {{\mathfrak {m}}}^{j})\) and \(\bar{\mathtt{a}}_i\big (\frac{\pmb {{\mathfrak {m}}}^j + \pmb {{\mathfrak {e}}}^j}{2}\big )\) using \(\{ +\pmb {\xi }^{j}(\omega _{k}) \}_{j=\ell }^{J}\) resp.  \(\{ -\pmb {\xi }^{j}(\omega _{k})\}_{j=\ell }^{J}\) with initial condition \(\pmb {{\mathfrak {m}}}^{\ell }=\pmb {g}\).

      2. (d)

         Compute \(H(t_{\ell },\pmb {g}_{\bar{h},i,l}^{+},\omega _{k})\) resp.  \(H(t_{\ell },\pmb {g}_{\bar{h},i,l}^{-},\omega _{k})\) in (5.4) based on \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{+,\ell ,k}\) resp.  \(\mathcal {S}_{\pmb {{\mathfrak {m}}}}^{-,\ell ,k}\) to determine \(d_{l,i}H(t_{\ell }, \pmb {g},\omega _{k})\) in (5.5) and store it.

   3. (6)

      Approximate \(\mathbb {E}[d_{l,i}H(t_{\ell },\pmb {g})]\) in (5.5) via Monte-Carlo estimation along with the variance reduction method of antithetic variates.

7. (VII)

   Set \(\nabla _{\mathcal {M}} w(t_{\ell },\pmb {g}) \approx \bigl ( d_{1}H(t_{\ell },\pmb {g}),\ldots ,d_{N}H(t_{\ell },\pmb {g}) \bigr )^{\top }\) with \(d_{i}H(t_{\ell },\pmb {g})\) as in (5.5), and compute \(\bar{\mathbf{u}}(t_{\ell },\mathbf{g}^\ell )\) as in (3.15), and \(\mathbf{m}^{\ell +1}\) via scheme (4.3b) with \(\mathrm{a}_{i}\bigl (\mathbf{g}^{\ell },\bar{\mathbf{u}}(t_{\ell },\mathbf{g}^{\ell })\bigr )\).

We now approximate the optimal pair \((\mathbf{m}^*, \mathbf{u}^*)\) by piecewise constant processes (indexed by the time interval [0, T]), whose value on \([t_j, t_{j+1}),\, j\in \{ 0,1,\ldots , J-1\}\) are respectively \(\mathbf{m}^{j}\) and \(\bar{\mathbf{u}}(t_{j},\mathbf{m}^j)\), where the iterates \(\{\mathbf{m}^{j}\}_{j=0}^{J} \) and \(\{\bar{\mathbf{u}}(t_{j},\mathbf{m}^j)\}_{j=0}^{J}\) are computed via Algorithm 5.1. To be more precise,

We computationally study the behavior of the optimal state and control for the switching dynamics of an ensemble of N particles by using the algorithms from Sect. 5. For this purpose, we employ discretely distributed random numbers from the GNU Scientific Library [10]. All computations are performed on an Intel Core i5-4670 3.40GHz processor with 16GB RAM in double precision arithmetic. The arising linear algebraic systems are solved by the Gaussian elimination method [10].

6.1 Test Studies

We start with test problems to compare the two methods from Sect. 5.2. For this purpose, we omit certain energy contributions in (1.1), and allow only one or two spins such that an exact solution of (3.12) becomes available.

Test problem 1: Consider the controlled problem for a single spin (\(\mathbf{J}=\mathbf{0}\)) of an isotropic material (\(\mathbf{D}=\mathbf{0}\)), and \(\delta =0\) in the cost functional; all other parameters \((C_\mathtt{ext}, \lambda , \nu , \alpha )\) are equal to 1. Then (3.12) is the backward heat equation

We shall use spherical harmonics to describe the exact solution of (6.1). Note that for any \(\pmb {m}= ({m}_1, {m}_2, {m}_3)\in \mathbb {S}^2\), the spherical harmonic \(w_{0,1}(\pmb {m})= {m}_3\) is an eigen-function of the Laplace-Beltrami operator with eigenvalue \(-2\) [9, Lemma 4.3.26], i.e.,

As \(w= \exp (-W)\), we know that w has to be positive, while \(w_{0,1}\) may also take negative values. Therefore we also use the constant spherical harmonic \(w_{0,0}(\pmb {m}):=1\). Consider the problem (6.1) with final time condition \(w(T,\pmb {m})=w_{0,1}(\pmb {m}) + 2\). We obtain this terminal condition by choosing the terminal payoff \(h(\pmb {m})= - \frac{1}{2}\log \big (w_{0,1}(\pmb {m}) + 2\big )\). Then the solution of (6.1) is

Moreover, we have the explicit formula for \(\nabla _{\mathbb {S}^2} w(t,\pmb {m})\), and hence for \(\nabla _{\mathbb {S}^2} W(t,\pmb {m})\):

Let \(\bar{\mathbf{u}}_\mathtt{exct}(t,\pmb {m})\) resp. \(\bar{\mathbf{u}}_\mathtt{app}(t,\pmb {m})\) be the function defined in (3.15) associated to (6.2) resp. (4.2), and \(\mathbf{m}_\mathtt{exct}^*(t)\) resp. \(\mathbf{m}_\mathtt{app}^*(t)\) be the solution of (2.3) with \(\bar{\mathbf{u}}_\mathtt{exct}\big (t,\mathbf{m}_\mathtt{exct}^*(t)\big )\) resp. \(\bar{\mathbf{u}}_\mathtt{app}\big (t,\mathbf{m}_\mathtt{app}^*(t)\big )\). By denoting the error

we show the behavior of \(\mathtt{err}(t)\) for different values of Monte-Carlo realizations with sample size M. For the simulation, we use \(T=0.5\), \(\tau =10^{-2}\), \(\bar{h}= \frac{1}{\sqrt{M}}\), \(\bar{\mathbf{m}}= \mathbf{e}_1\), and other parameters as specified at the beginning of this subsection. We observe that the error \(\mathtt{err}(t)\) for Method B is significantly smaller (by a factor of \(\frac{1}{20}\) in our simulations) if compared to Method A, see Fig. 1. Moreover, at least \(M\approx 10^6\) realizations are needed to balance the approximate computation via Method B with the remaining error sources.

Remark 6.1

Computational studies with respect to the parameters \(\big (\tau ,\bar{h},M\big )\) show that, independent of \(\tau \), it is beneficial to choose \(\bar{h} = \mathcal {O}\big (\frac{1}{\sqrt{M}}\big )\) to approximate the \(\nabla _\mathcal {M} w\) (hence \(\nabla _\mathcal {M} W\)) accurately. For choice \(\bar{h} \ll \mathcal {O}(\frac{1}{\sqrt{M}})\), irrespective of the Method A or Method B, we observe a strongly oscillatory behavior of the solution when \(M\le 10^4\).

Test problem 2: We study the interaction of two isotropic (\(\mathbf{D}=\mathbf{0}\)) spins for \(\alpha =0=\delta \), and other parameters \((C_\mathtt{ext}, \lambda , \nu )\) are equal to 1. Let us first recall how the spherical harmonics on a single sphere \(\mathbb {S}^2\) generalize to the manifold \(\mathcal {M}=(\mathbb {S}^2)^N= (\mathbb {S}^2)^2\) most naturally. Indeed, because \(\mathcal {M}\) is a tensor product of spheres, and the spherical harmonics form an orthogonal basis on the single sphere, the tensor products of spherical harmonics form an orthogonal basis on \(\mathcal {M}\). It is therefore reasonable to expect that the simplest meaningful test problems on \((\mathbb {S}^2)^2\) can be constructed with terminal time conditions which are products of low-order spherical harmonics and which are eigen-functions of the transformed Bellman equation. Because of the spin interaction, the first order coefficient \(\mathbf{b}\) in (3.12) does not vanish any more. Therefore, we combine the functions \(w_{0,0}(\pmb {m}_1) w_{0,0}(\pmb {m}_2), w_{0,0}(\pmb {m}_1)w_{0,1}(\pmb {m}_2)\) and \(w_{0,1}(\pmb {m}_1)w_{0,0}(\pmb {m}_2)\) with \(\pmb {m}_1=({m}_{1,1},\, {m}_{1,2},\,{m}_{1,3})\in \mathbb {S}^2\) and \(\pmb { m}_2=({m}_{2,1},\, {m}_{2,2},\,{m}_{2,3})\in \mathbb {S}^2\) to pose a test problem on \((\mathbb {S}^2)^2\). Denoting by

we consider the following version of (3.12),

with the positive semi-definite matrix \(\mathbf{J}\): for any \(\mu >0\),

Since \(\alpha =0\), we have

We compute the tangential gradient of the functions \(w_{00,00},\,w_{01,00}\) and \(w_{00,01}\):

Observe that

and \(w_{01,00}\) and \(w_{00,01}\) are eigen-functions of the Laplace-Beltrami operator on \((\mathbb {S}^2)^2\) with eigenvalue \(-2\). Thus the exact solution of (6.3) is given by

Moreover, we compute \(\nabla _{(\mathbb {S}^2)^2} W\big (t, \pmb {m}_1, \pmb {m}_2\big )\) as

Note that the choice of W corresponds to the terminal payoff

Similar to test problem 1, we define \(\mathtt{err}(t)\) and study its behavior in time t for different Monte-Carlo realizations with sample size M by using Method B, see Fig. 1c. The simulation is made for the following choice of parameters: \(T=0.5\), \(\tau =10^{-2}\), \(\bar{h}= \frac{1}{\sqrt{M}}\), \(\bar{\mathbf{m}}=\big (\mathbf{e}_1, \mathbf{e}_2\big )\), and other parameters as specified in the problem. We observe that the error \(\mathtt{err}(t)\) decreases if one increases the sample size M.

Time evolution of a single trajectory of the error \(t\mapsto \mathtt{err}(t)\) : a by using Method A, b by using Method B for test problem 1, and c by using Method B for test problem 2. Note the different ranges on the axes

Full size image

We observe that the error \(\mathtt{err}(t)\) for the both test problems 1 and 2 is of the same magnitude as the error made in the approximation of \(\nabla _\mathcal {M} w\) (hence \(\nabla _\mathcal {M} W\)).

Optimal control of two interacting isotropic spins. Remaining in the setting of test problem 2 we next study the time evolution of a single trajectory of the optimal state, as well as the magnitude and direction of the optimal control. In this case, the trajectory of the optimal control lies in \(x_1\)-\(x_2\) plane to balance the random influences; see Fig. 2.

Test problem 2: time evolution of a single trajectory of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), and the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\) for \(i=1,2\) (Color figure online)

Full size image

Remark 6.2

Computational studies for both test problems suggest stability of the scheme (4.3a)–(4.3b). However, convergence resp. termination of the scheme depends crucially on the given parameters in Problem 1.1. For choices

an exponential overflow occurs during truncation in simulations, and therefore the computed value of \(w\big (t,\pmb {m}\big )\) in (4.2) is set to zero then. Hence, in this case, \(\log \big (w(t,\pmb {m})\big )\) is not defined, and thus the approximation procedure to approximate \(\nabla _{\mathcal {M}} W(t,\pmb {m})\) terminates. This is one reason that Examples 5.1 and 5.2 from [7] may not directly be simulated here. Notice that no exponential overflow occurs for both test problems above, since \(\delta = 0\).

6.2 Optimal Control of Three Interacting Spins

We now study an ensemble of \(N=3\) particles, which additionally are subjected to exchange forces. We are interested in the switching control for one (\(i=2\)) of these particles from \(\bar{\mathbf{m}}_{2}\) (at initial time) to \(-\bar{\mathbf{m}}_{2}\) at given final time T. Take \(h(\mathbf{m}(T))= \frac{1}{2}\Vert \mathbf{m}(T)- \widetilde{\mathbf{m}}(T)\Vert _{(\mathbb {R}^3)^N}^2\), where the deterministic target profile \(\widetilde{\mathbf{m}}: [0,T]\rightarrow (\mathbb {S}^2)^3\) is given by

We use again Method B to approximate \(\nabla _\mathcal {M} W\). To simulate the optimal pair of the underlying problem, we have used the methodology described in Sects. 4.1 and 5.1–5.3, along with the following set of parameters:

with the positive semi-definite matrix \(\mathbf{J}\) such that for any \(\mathbf{m}=(m_1, m_2,\ldots , m_N)\in (\mathbb {S}^2)^N\)

Time evolution of a single trajectory of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\), and the angle between optimal pair \(t\mapsto \ll u_i^*(t), m_i^*(t)\gg \) for \(i=1,2,3\) (Color figure online)

Full size image

In this case, the minimum value of the cost functional is \(\mathcal {J}_\mathtt{sto}^*\equiv 0.9078\). Though the first and third spins start already at the desired state, it is due to the noise, and the exchange forces in particular, that the optimal control is acting on the whole time interval and on all spins. For the second spin, we observe that at the beginning and end, less control is needed opposed to the applied control at the intermediate times; see Fig. 3. The orthogonality of the optimal pair \((\mathbf{m}^*, \mathbf{u}^*)\) (e.g. Remark 3.4) is shown in Fig. 3g–i by displaying the temporal evolution

6.3 Optimal Control of Four Interacting Spins

We consider here the switching control for an ensemble of \(N=4\) particles.

Set-up 1: We use the parameters as in Sect. 6.2 with \( \bar{\mathbf{m}}=\big ( \mathbf{e}_1, -\mathbf{e}_1, \mathbf{e}_{\mathbf{1}}, -\mathbf{e}_1\big )\), and \( \widetilde{\mathbf{m}}(t)=\big (\mathbf{e}_1,\widetilde{\mathbf{m}}_2(t), \mathbf{e}_1, \widetilde{\mathbf{m}}_2(t)\big )\). In this case, the first and third spins start already at the desired state; the associated optimal controls are acting on the whole time interval. Moreover, for the second and fourth spins, significant controls are required to approximately meet the terminal state profile. The time evolution of \(t\mapsto \Vert u_2^*(t)\Vert _{\mathbb {R}^3}^2\) is similar to the results for \(N=3\) spin constellations (see Fig. 3e), while \(\Vert u_4^*(t)\Vert _{\mathbb {R}^3}^2\) is delayed in time for the fourth spin. We observe a loop of the orientation of \(u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}~(i=2,4)\) close to the terminal time; see Fig. 4.

Time evolution of a single trajectory of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), and the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\) with set-up 1 for \(i=1,2,3,4\) (Color figure online)

Full size image

Set-up 2: We use same parameters as in set-up 1 with \( \bar{\mathbf{m}}=\big ( \mathbf{e}_1, -\mathbf{e}_1, -\mathbf{e_1}, \mathbf{e}_1\big )\) and \(\widetilde{\mathbf{m}}(t)=\big (\mathbf{e}_1,\widetilde{\mathbf{m}}_2(t), \widetilde{\mathbf{m}}_2(t), \mathbf{e}_1\big )\). For the second and third spins, significantly synchronous controls at intermediate times are required to meet approximately the desired target profile. Like in Set-up 1, we also observe the formation of loops of the orientation of \(u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}~(i=2,3)\) close to the terminal time; see Fig. 5.

Time evolution of a single trajectory of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), and the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\) with set-up 2 for \(i=1,2,3, 4\) (Color figure online)

Full size image

6.4 Optimal Control of Ten Interacting Spins

We consider here an ensemble of \(N=10\) particles to optimally control the dynamics to reach a deterministic target profile

within finite time T at minimized expected external energy with initial configuration

To simulate the optimal pair of the underlying problem, we take again \(\mathbf{D}=\mathbf{0}\), \(\mathbf{J}\) as in (6.5), \(h(\mathbf{m}(T))= \displaystyle \frac{1}{2} \Vert \mathbf{m}(T)- \widetilde{\mathbf{m}}(T)\Vert _{(\mathbb {R}^3)^N}^2\), and the following set of parameters:

In Fig. 6a, we visualize the behavior of the optimal state \(\mathbf{m}^*\). Due to the large damping coefficient \(\alpha =1.0\), we observe fast switching dynamics of the optimal state. With the choice \(\lambda =1\) the control is penalized more strongly than in the previous experiments, which has a noticeable effect on the magnitude of \(\mathbf{u}^*\), compare Figs. 5e–h and 6c, e. At the beginning a stronger control is applied to move towards the desired target profile. Because of the large noise intensity \(\nu \), and the less control, some particles of this ensemble do not reach the target profile appropriately. For illustration, we plotted the behavior of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), and the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\) for \(i=3,7\); see Fig. 6.

A single realization of the optimal state \(m^*\) (red) at final time T, the time evolution of the optimal state \(t\mapsto m^*_i(t)\) (red), the direction of the optimal control \(t\mapsto u_i^*(t)\Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{-1}\) (blue), and the magnitude of the optimal control \(t\mapsto \Vert u_i^*(t)\Vert _{\mathbb {R}^3}^{2}\) for \(i=3,7\) (Color figure online)

Full size image

Let us specify the cost of the proposed method in terms of computations and storage. In the definition of \({{\mathcal {J}}}_\mathtt{sto}\) resp. w in Problem 1.1 resp. (4.2), we approximate the expectation via Monte-Carlo estimation along with the variance reduction method of antithetic variates. There, each realization of the integrals requires to store the iterates \(\{ \mathbf{m}^j\}_{j=0}^J\) resp. \(\{\pmb {{\mathfrak {m}}}^j \}_{j=0}^J\) with \(\mathcal {O}(NJ)\) storage complexity due to the necessity of piecewise interpolation as in (5.1). The latter iterates are computed via scheme (4.3a)–(4.3b) by solving 2N linear \(3\times 3\) systems, which can be done analytically in \(\mathcal {O}(1)\) time and in parallel, cf. the discussion in [17, Section 2.3]. The most computationally intensive part is to approximate \( \nabla _\mathcal {M} w(t_{j},\mathbf{m}^j)\), whose accuracy has a direct impact on both, the optimal solution and control. Comparative computational studies show that at least \(\mathcal {O}({10}^6)\) realizations are needed to approximate \(\nabla _\mathcal {M} w(t_{j},\mathbf{m}^j)\), whereas \(\mathcal {O}({10}^4)\) realizations are sufficient to approximate the cost function. This dependence on the sample size \(M\in \mathbb {N}\) to approximate the optimal control is strengthened in the case of fast changing optimal states.