Pontryagin Neural Networks for the Class of Optimal Control Problems With Integral Quadratic Cost

Introduction

Optimal control theory has been a major area of interest in various scientific fields throughout the last century. Optimal Control Problems (OCPs) are fundamental in many disciplines, leading to the development of numerous numerical methods aimed at optimising cost functions that are central to OCPs. In general, OCPs can be addressed using two main approaches: direct methods and indirect methods [1].

Direct methods involve discretising the continuous states and controls to transform the continuous problem into a Nonlinear Programming (NLP) problem [2–4]. This transformation results in a finite constrained optimisation problem solvable by various algorithms designed to find local minima, such as trust-region methods [5]. Direct methods have been applied to a wide range of OCPs, particularly in the aerospace sector (e.g., [6–9]). However, the general NLP problem is considered NP-hard, meaning that the computational time required to find the optimal solution cannot be predetermined and can be substantial. This uncertainty in convergence time can undermine the reliability of these methods, especially in real-time applications.

To address these challenges, researchers have explored transforming OCPs from non-convex formulations into convex optimisation problems [10, 11]. Convex problems are computationally tractable because their associated numerical algorithms guarantee convergence to a global optimum in polynomial time. This approach involves applying convexification techniques to reformulate the optimal guidance problem into a convex one. Such methodologies have been successfully used in various aerospace applications, including planetary landing [10, 11], atmospheric entry guidance [12, 13], rocket ascent guidance [14], and low-thrust trajectory optimisation [15].

Indirect methods, on the other hand, derive first-order necessary conditions by applying the Pontryagin Minimum Principle (PMP) or the calculus of variations. These necessary conditions lead to a Two-Point Boundary Value Problem (TPBVP) that involves state and costate variables, which must be solved numerically. Available numerical methods for indirect approaches include single and multiple shooting methods [16, 17], orthogonal collocation [18], and pseudospectral methods [19]. The commercial software GPOPS-II, based on the work presented in [20], represents the current state-of-the-art. While the indirect method guarantees an optimal solution, solving the TPBVP can be challenging due to the sensitivity of the solution to initial guesses, particularly for costate variables that lack physical interpretation. Although direct methods may not always yield optimal solutions, they are generally easier to compute and are widely used in the scientific community.

In principle, these two approaches can be combined. For instance, a solution obtained via a direct method can serve as an initial guess for an indirect method, potentially refining the solution and ensuring optimality.

In this work, we addressed a class of OCPs with integral quadratic cost functions by applying the PMP and solving the resulting TPBVP using a specialised Physics-Informed Neural Network (PINN) framework known as the Extreme Theory of Functional Connections (X-TFC) [21–24]. In our approach, PINNs are tailored and trained to learn optimal control actions by enforcing the physics constraints represented by the TPBVP arising from the application of the PMP. We refer to these specialised PINNs as Pontryagin Neural Networks (PoNNs).

The main contributions and objectives of this paper are threefold: 1) to demonstrate the feasibility of using PINN-based frameworks to directly learn TPBVP solutions and, consequently, open-loop optimal control with high accuracy and low computational time; 2) to provide an estimate of the generalisation error associated with the proposed methodology; and 3) to compare the performance of PoNNs with the commercial software GPOPS-II [20].

Physics-Informed Neural Networks (PINNs) are a machine learning framework that incorporates physical laws, expressed as differential equations (DEs), into the training of neural networks (NNs) by including them as regularisation terms in the loss function [25]. In classical PINNs, the DE constraints are included in the loss function alongside the residuals computed within the domain. However, this approach has a significant drawback: it requires the simultaneous satisfaction of the DEs and their associated boundary conditions (BCs) or initial conditions (ICs) during training, which can complicate the optimisation process [26].

The X-TFC framework enhances the classical PINN approach by using Constrained Expressions (CEs) from the Theory of Functional Connections (TFC) [27]. The general form of a CE is given by

where $f(\mathbf{x})$ is the unknown function to be approximated, $f_{CE}(\mathbf{x},g(\mathbf{x}))$ is the CE, $A(\mathbf{x})$ is a functional that analytically satisfies the given linear constraints, and $B(\mathbf{x},g(\mathbf{x}))$ projects the free function $g(\mathbf{x})$—which exists over the constraints—onto the space of functions that vanish at the constraints [28].

TFC is a functional interpolation technique applicable to various mathematical problems, including the solution of differential equations (DEs) [29] and the addressing of OCPs via the PMP approach. While the original TFC employs linear combinations of orthogonal polynomials as the free function, this can lead to computational challenges, especially for large-scale partial differential equations (PDEs), due to the curse of dimensionality [29]. To mitigate these issues, PINN TFC-based methods use neural networks as the free function. Specifically, X-TFC employs Extreme Learning Machines (ELMs), which are shallow neural networks with random features, also known as Random Projection Networks.

X-TFC has been successfully applied to solve integro-differential equations [30] with applications in rarefied gas dynamics [31, 32] and radiative transfer equations [33]. It has also been used for stiff systems of ODEs in nuclear reactor dynamics [34] and stiff chemical kinetics [35], as well as inverse problems for system identification [36, 37] and parameter estimation in epidemiological models [38]. Additionally, X-TFC has been applied to a class of OCPs for aerospace applications using both indirect methods [39, 40] and the Bellman Optimality Principle [41].

This paper is organised as follows. First, we describe how PoNNs are constructed and trained to learn control actions for the class of OCPs with integral quadratic cost functions. We then provide an upper bound estimate on the generalisation error of PoNNs in learning solutions for the considered OCP class, in terms of training error and the number and selection of training points. Two benchmark OCPs—one linear and one nonlinear—are used to demonstrate the effectiveness of PoNNs in learning optimal control actions and to validate our theoretical findings. We also compare the performance of PoNNs with the commercial software GPOPS-II. Finally, we present our conclusions.

Materials and Methods

Pontryagin Neural Networks (PoNNs) are a specialised class of Physics-Informed Neural Networks (PINNs) designed to solve Optimal Control Problems (OCPs) through the application of the Pontryagin Minimum Principle (PMP). In this study, we focused on OCPs characterised by integral quadratic cost functions.

The term Pontryagin Neural Networks (PoNNs) was introduced to succinctly describe PINNs that are tailored to learn optimal control strategies via PMP [39]. The unknown solutions of the resulting Two-Point Boundary Value Problem (TPBVP), specifically the states and costates, are represented using a PINN framework based on the Extreme Theory of Functional Connections (X-TFC), developed by the authors.

In PINN TFC-based frameworks, unknown solutions of differential equations are modelled using Constrained Expressions (CEs) from the Theory of Functional Connections (TFC), with neural networks serving as free functions. When Deep Neural Networks (DNNs) are employed as free functions, the framework is referred to as Deep-TFC [28], where gradient-based optimisation methods like the ADAM optimiser [42] or the L-BFGS optimiser [43] are used to adjust the parameters.

Alternatively, when shallow neural networks trained using the Extreme Learning Machine (ELM) algorithm [44] are used as free functions, the framework is known as Extreme-TFC (X-TFC) [21]. The ELM algorithm initialises the input weights and biases randomly, keeping them fixed during training, which means that only the output weights are optimised [44]. This approach allows for robust and fast training using least-squares methods.

For linear differential equations, X-TFC optimises the free-function parameters (i.e., the neural network output weights) using a single-pass least-squares method. In the case of nonlinear differential equations, iterative least-squares procedures are required to optimise the parameters, as detailed in [21].

This section outlines the design and training process of PoNNs to learn optimal control actions for the class of OCPs under consideration. Additionally, following Refs. [22, 45], we derive an estimate for the upper bounds on the generalisation error of PoNNs in learning solutions for this class of OCPs.

A general OCP can be formulated as:

Here, $\mathit{x}\in {\mathrm{\Omega }}_x\subseteq {\mathbb{R}}^n$ represents the states that may be subject to path constraints, $\mathit{u}\in \mathcal{U}\subseteq {\mathbb{R}}^m$ denotes the controls that may be subject to inequality constraints, and $\mathcal{C}\subset {\mathrm{\Omega }}_x$ defines the terminal conditions. The cost function includes both a Meyer term (dependent on the final state and time) and a Lagrangian term (dependent on the states and controls over time). The final time $t_f$ can be fixed or a function of the initial states ${\mathit{x}}_0$ (e.g., in time-free problems). Applying PMP transforms the OCP into the following TPBVP:

with $t\in [t_0,t_f]$, where $H$ is the Hamiltonian of the problem, and $H_{\mathit{u}}=\frac{\partial H}{\partial \mathit{u}}$, $H_{\mathit{x}}=\frac{\partial H}{\partial \mathit{x}}$, $H_{\mathit{\lambda }}=\frac{\partial H}{\partial \mathit{\lambda }}$, with $\mathit{\lambda }$ being the costates. Equation 2 encompasses all possible transversality conditions for unconstrained OCPs, although the specific conditions applied depend on the problem at hand. Further details are available in [46, 47].

In this work, we focus on unconstrained OCPs with a fixed final time and a fixed final state, leading to the following TPBVP:

Here, $\mathit{x}$, $\mathit{\lambda }$, $H_{\mathit{\lambda }}$, $H_{\mathit{x}}$, ${\mathit{x}}_0$, and ${\mathit{x}}_f$ are all $n$-dimensional vectors, with $n$ being the number of states. Equation 3 is well-posed and can be expressed compactly as:

where $\mathit{X}={[{\mathit{x}}^{\top },{\mathit{\lambda }}^{\top }]}^{\top }\in {\mathbb{R}}^{2n}$ and $\mathit{F}={[H_{\mathit{\lambda }}^{\top },-H_{\mathit{x}}^{\top }]}^{\top }\in {\mathbb{R}}^{2n}$. We assume that $\mathit{F}$ is globally Lipschitz continuous, meaning that there exists a constant $C_f\ge 0$ (independent of $t$) such that:

To apply the X-TFC algorithm to solve Equation 4, we need to specify the training set $\mathcal{S}$ and define the residuals. Considering the time domain $[t_0,t_f]\subset \mathbb{R}$, we selected the training set $\mathcal{S}\subset [t_0,t_f]$ based on suitable quadrature points, such as those corresponding to a composite Gauss quadrature rule or randomly selected points, as suggested in Ref. [45]. The training points are defined as:

Using the constrained expressions, we approximate the solution of Equation 4. The states are approximated by:

The costates are approximated by:

where

and

In these expressions, $\mathit{w}={[w_1,\dots ,w_L]}^{\top }$ and $\mathit{b}={[b_1,\dots ,b_L]}^{\top }$ are the input weights and biases, respectively, and $L$ is the number of hidden neurons. Since $g_{x_i}$ and $g_{{\lambda }_i}$ are single-layer neural networks with random features, the weights $\mathit{w}$ and biases $\mathit{b}$ are randomly initialised and remain fixed during training; only the output weights $\mathit{\beta }$ are optimised. Notably, $x_i(t;{\mathit{\beta }}_{x_i})$ and ${\lambda }_i(t;{\mathit{\beta }}_{{\lambda }_i})$ constitute the PoNNs.

We define the residuals as:

where

with

To train the PoNNs and learn the unknown solutions, we minimise the following vector-valued loss function:

With the loss function defined, the X-TFC algorithm is fully specified and can be executed to learn the solution of the TPBVP in Equation 4 resulting from the application of the PMP. The learnt PoNNs are denoted as ${\mathit{X}}^{*}={\mathit{X}}_{{\mathbf{\Xi }}^{*}}$, where ${\mathbf{\Xi }}^{*}$ are the optimised output weights obtained from training.

If the loss function in Equation 7 is linear, it simplifies to a system of linear algebraic equations of the form $H\mathbf{\Xi }=\mathit{T}$, which can be solved using a single-pass least-squares method. For nonlinear loss functions, an iterative least-squares procedure is necessary, as described in [21].

Figure 1 illustrates the flowchart for learning optimal control problem solutions using the proposed PoNNs.

Figure 1. Flowchart for learning optimal control problem solutions using the proposed PoNNs.

Estimation of the Generalisation Error

We are interested in estimating the generalisation error for X-TFC, which is defined as

where ${\epsilon }_G$, according to [45] and references therein, is the error emerging from approximating the true solution, $\mathit{X}(t)$, via X-TFC, within the domain $t_0\le t\le t_f$. We estimate the generalisation error in terms of the training error defined as,

which is computed a posteriori once the training is completed.

According to Ref. [45], we also need the following assumptions on the quadrature error: for any function $\mathit{y}\in C^l([t_0,t_f])$, the quadrature rule corresponding to the quadrature weights $w_i$ at points $t_i\in \mathcal{S}$, with $1\le i\le N$, satisfies,

where $\alpha >0$ and different order quadrature rules can be used (details can be found in [45]). The estimate of the generalisation error for X-TFC is given by the following theorem.

Theorem 1. Let $\mathit{X}\in C^k([t_0,t_f])$ be the unique classical solution of the general TPBVP (Equation 4), resulting from the PMP application to Equation 1, with the dynamics $\mathit{F}$ satisfying Equation 5. Let ${\mathit{X}}^{*}={\mathit{X}}_{{\mathbf{Xi}}^{*}}$ be the solution approximated via X-TFC, corresponding to the loss function (Equation 7). Then, the generalisation error (Equation 8) can be estimated as,

where ${\epsilon }_T$ is given by Equation 9, the constant $C_{DE}$ is given by,

and

which is a function of the residuals $\mathcal{R}$,

Also

which is a positive constant depending on the number of training points, the quadrature scheme used, and the residuals evaluated on the training points [45].

Proof

Using the definitions of the residuals (Equation 6), and the system (Equation 4), we can verify that the error $\widehat{\mathit{X}}={\mathit{X}}^{*}-\mathit{X}$ satisfies the following equation,

Here, we have denoted $\mathcal{R}={\mathcal{R}}_{{\mathbf{\Xi }}^{*}}$ for convenience of notation. Multiplying both sides of the system (Equation 11) by $\widehat{\mathit{X}}$ yields

where $\cdot$ is the inner product. That is,

By leveraging on the following inequality ${\left(\widehat{\mathit{X}}-\mathcal{R}\right)}^2\ge 0$, Equation 12 is bounded by,

By Equation 5, we have the following,

Rearranging we get,

Integrating the above inequality over $[t_0,\bar{T}]$ for any $t_0\le \bar{T}\le t_f$, we obtain,

where,

with $\Vert \widehat{\mathit{x}}(t_0){\Vert }^2=0$ because of the CE.

Also,

since $\Vert \mathcal{R}{\Vert }^2\ge 0$ and $t_0\le \bar{T}\le t_f$. Thus,

Applying the integral form of the Grönwall’s inequality to the above, we get

Integrating over $[t_0,t_f]$ in $d\bar{T}$,

That is,

Thus,

Finally,

where,

This completes the proof.

The estimate (Equation 10) indicates that the generalisation error when training PoNNs is small under several conditions. The training error for the PoNN must be reasonably small (e.g., ${\epsilon }_T\ll 1$). Moreover, the training error can only be computed a posteriori, as we have no prior control over it. The error associated with the training points depends on the choice and the number of training points $N$ as well as on the training constant $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}$. The latter also depends on the residual of the PoNN, ${\mathit{X}}^{*}$, and thus, indirectly, on the number of training points $N$. That is, $N$ needs to be large enough so that $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}^{\frac{1}{2}}{N}^{-\frac{\alpha }{2}}\ll 1$. The constant $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}$ also depends on the architecture of the underlying PoNN. The $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}$ evaluation depends on the features of the governing DEs, and the collocation of the training points cannot be worked out in the setting of Theorem 1 [48]. The constant $C_{\mathrm{p}\mathrm{d}\mathrm{e}}$ encodes the stability of the DEs and depends on both the exact (unknown) solution $\mathit{X}$ and the trained PoNN, ${\mathit{X}}^{*}$, that needs to be bounded. For OCPs with fixed final state, the term $\widehat{\mathit{\lambda }}(t_0)$ can be computed for linear problems. In this case it is related to the transition matrix. For nonlinear problems, $\widehat{\mathit{\lambda }}(t_0)$ cannot be known. If it were known, the main limitation of the indirect method for solving OPCs would be removed (e.g., the sensitivity of the guess on the costates initial value). Moreover, if this were the case, we would solve the problem with a perfect shooting method. However, if the domain $[t_0,t_f]$ is compact, we can guarantee that if the other terms in the ${\epsilon }_G$ are small enough, then $\widehat{\mathit{\lambda }}(t_0)$ will also be small enough (e.g., $\widehat{\mathit{\lambda }}(t_0)\ll 1$). Therefore, PoNN ${\mathit{X}}^{*}$ will generalise well the unknown solution $\mathit{X}$ (e.g., ${\epsilon }_G\ll 1$). Conversely, if the OCPs have a free final state, in the TPBVP (Equation 4) we will have the condition $\mathit{\lambda }(t_f)={\mathit{\lambda }}_{\mathit{f}}$ instead of $\mathit{x}(t_f)={\mathit{x}}_{\mathit{f}}$. In this case, the condition on the costate at the final time can be considered as an initial condition for a backward problem, and thus the TPBVP is tackled as a first-order Initial Value Backward Problem.

The main point of the estimate (Equation 10) is to provide an upper bound of the generalisation error in terms of the training and quadrature errors, regardless of the choice of the PoNN architecture. The estimate (Equation 10) holds for any expansion of $\mathit{X}$, including NNs (of any type). There is no guarantee that the training error will be small if NNs are used to approximate $\mathit{X}$. However, it can be expected that if the DE is stable and the training points collocation scheme is accurate (e.g., there is some control over $C_{\mathrm{p}\mathrm{d}\mathrm{e}}$ and $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}$), the training error will be small, leading to a small generalisation error. Nevertheless, for some PoNN architectures, where the training error was small, but $C_{\mathrm{p}\mathrm{d}\mathrm{e}}$ and $C_{\mathrm{q}\mathrm{u}\mathrm{a}\mathrm{d}}$ were not (e.g., for high $t_f$ and/or not enough training points $N$), the overall estimate would not be small enough (e.g., ${\epsilon }_G$ would not be $\ll 1$). In this scenario, PoNN ${\mathit{X}}^{*}$ would not generalise well the unknown solution $\mathit{X}$.

Results

To test our theoretical findings, two benchmark OCPs are studied: 1) linear hypersensitive Feldbaum problem and 2) nonlinear hypersensitive Feldbaum problem.

All the selected problems are coded in MATLAB 2022a and run on an Intel Core i7 - 9700 CPU PC with 64 GB of RAM.

Problem 1: Linear Hypersensitive Feldbaum Problem

The problem considered here is a linear hypersensitive OCP taken from [49, 50], also known as the Hypersensitive Feldbaum Problem (HFP). The problem is cast as follows:

The resulting TPBVP after the PMP application is,

The analytical solution is available and provided in [50].

For the PoNN algorithm, we will use the following hyperparameters in all cases: Gaussian activation function, input and bias are sampled from $\mathcal{U}(-3;3)$. It should also be noted that the computed times of execution are reported on their own scale.

In Figure 2, the two plots on the left show the analytical, PoNNs, and GPOPS-II solutions. On the plots on the right, the absolute errors of PoNNs and GPOPS-II at the points when compared with the analytical solution are shown. In the first row for the State variable, in the second row for the Control one. We used $L=100$, which was trained on 100 randomly generated points. GROPS-II used $N=61$ points.

Figure 2. Computed solutions and absolute errors for the linear test problem.

Figure 3 shows the comparison with the results obtained by GPOPS-II. We report convergence results when increasing the number of neurons $L$ fixing $N=61$ and using exactly the points given by the (adaptive) procedure of GROPS-II. We note that in the result, the error made by the commercial software is strongly affected by an anomalous value obtained at the final time, see Figure 2.

Figure 3. Convergence of the PoNN algorithm when increasing the number of neurons $L$, comparison with the result obtained by GPOPS-II. The training points are fixed and coincide with those used by GPOPS-II $N=61$.

In Figure 4 we report the obtained results by applying our algorithm. It should be noted that the computation time can be affected by the (machine precision) lower rank.

Figure 4. Convergence history for the application of the PoNN method to the linear problem. At the top, is the case with $N=100$ fixed training points, at the bottom, the one with $L=100$ fixed number of neurons.

The numerical tests confirm the theoretical findings. The comparison against GPOPS-II shows that PoNNs outperform GPOPS-II both in terms of accuracy and computational time. It should be noted that GPOPS-II was run with the default options, i.e., we did not require machine-level accuracy. In all cases, PoNNs achieve the machine-level accuracy with a lower computational time than GPOPS-II. Moreover, this last only reaches a default accuracy. Therefore, if GPOPS-II were pushed to machine level accuracy the gap in computation time would be higher. This assesses the superiority of PoNNs in terms of both accuracy and computation time.

Problem 2: Nonlinear Hypersensitive Feldbaum Problem

The problem considered here is a nonlinear version of the HFP considered before. The problem is formulated as follows:

The resulting TPBVP after the PMP application is,

The analytical solution is not available. Thus, to compute the error, we estimate it over the residual on the training points. To validate the accuracy of the target PoNN solution, we have compared it with the one generated with GPOPS-II and also evaluated its residuals. In all simulations, the hyperparameters used are: Gaussian activation function, input and bias sampled from $\mathcal{U}(-10;10)$.

In Figure 5 we report state and control time evaluations on training points. Training is performed on $N=100$ randomly sampled points with $L=100$ hidden neurons for the PoNN method, while GPOPS uses $N=86$ training points.

Figure 5. Visual comparison of the numerical solution. PoNN on N = 100 random training points. GPOPS uses 86 points for training.

In Figure 6 we show the residuals obtained when the parameters are increased: at the top, the number of neurons is increased while the number of training points (randomly generated) is fixed at $N=100$; at the bottom the number of training points is increased while the number of neurons is fixed at $L=100$.

Figure 6. Convergence and stability results. On the ordinate is the computed residual, which is taken as an estimate of the overall error.

Even for the nonlinear case, the numerical tests confirm the theoretical findings.

Discussion

In this paper, we have introduced PoNNs and used them to learn the solutions to OCPs with integral quadratic cost, which were tackled via the indirect method.

The proposed approach can be applied to all research fields where OCPs with integral quadratic cost must be solved. The low computational time makes PoNNs attractive for real-time applications.

We have derived upper bounds on the generalisation error, in terms of the training error and number and choice of training points, of PoNNs in learning OCP solutions. The theoretical results are validated with numerical experiments on a benchmark linear and a benchmark nonlinear OCPs. Furthermore, the results are compared with those generated with the commercial software GPOPS-II. PoNNs outperform the GPOPS-II both in terms of accuracy and computation time. It should be noted that the GPOPS-II was run with the default options, i.e., we did not require machine-level accuracy. In all cases, PoNNs reach a higher level of accuracy with a lower computation time than GPOPS-II reaches the default accuracy. Therefore, if GPOPS-II were pushed to machine level accuracy the gap in computation time would be higher. This assesses the superiority of PoNNs in terms of both accuracy and computation time.

It should also be noted that PoNNs are not (at all) optimised as GPOPS-II is. For instance for a more complex problem, where a domain decomposition may be required to achieve better performance, GPOPS-II will most likely achieve better results than PoNNs. Therefore, leveraging on this paper’s results (in particular the superiority of PoNNs on simple benchmark problems such as those considered in this manuscript), work is in progress to optimise PoNNs. In particular, the authors are working on adaptive domain-decomposition techniques.

A different and complementary approach to solving OCP comes from the application of the Bellman Principle of Optimality (BPO). The latter provides the necessary and sufficient conditions for optimality which yields in writing the so-called Hamilton-Jacobi-Bellman Equation (HJB) equation. The HJB is a nonlinear PDE that plays a crucial role in the optimal control theory. The solution of the HJB equation is the value function in all the time-state space. The value function represents the optimal cost-to go for a dynamical system with an associated cost function. Once the value function is known, the optimal control is also known as it is the gradient of the value function to the state variables. This means that, if we learn the solution of the HJB equation in the entire time-state space, we would also consequentially learn the optimal control actions in the entire time-state space in a closed-loop fashion. However, since it is not trivial and, in many cases, it is most likely impossible to find an analytical solution to the HJB equation, one has to resort to numerical methods. The main issue in solving HJB PDEs is that the spatial complexity of the problem increases exponentially with respect to the number of dimensions of the system. Thus the majority of the state-of-the-art numerical methods to tackle PDEs struggle to solve the HJB equation. The reason why the HJB equation should be addressed in the future is that, unlike the PMP, which provides just a necessary condition for the optimality and results in an open-loop locally optimal trajectory, the HJB equation can guarantee a necessary and sufficient condition for the optimality and provides a closed-loop solution [51, 52].

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

ES and FC contributed to conception and design of the study. ES and DD wrote the first draft of the manuscript. All authors contributed to the article and approved the submitted version.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. FC was partially supported by PRIN 2022 PNRR P2022WC2ZZ “A multidisciplinary approach to evaluate ecosystems resilience under climate change” and by PRIN 2022 2022N3ZNAX “Numerical Optimization with Adaptive Accuracy and Applications to Machine Learning.”

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

FC and DD are members of the INdAM research group GNCS and were partially supported by GNCS-INdAM; this support is gratefully acknowledged.

References