Year: 2021

Abstract:
Synthesizing optimal controllers for dynamical systems often involves solving optimization problems with hard real–time constraints. These constraints restrict the class of numerical methods that can be applied: computationally expensive but accurate numerical routines often have to be replaced with fast and inaccurate methods, trading inference time for solution accuracy. This paper provides techniques to improve the quality of optimized control policies given a fixed computational budget. We achieve the above via a hypersolvers (Poli et al., 2020a) approach, which hybridizes a differential equation solver and a neural network. The performance is evaluated in direct and receding–horizon optimal control tasks, where the proposed approach shows consistent Pareto improvements in solution accuracy and control performance.

text: Synthesizing optimal controllers for dynamical systems often involves solving optimization problems with hard real–time constraints. These constraints restrict the class of numerical methods that can be applied: computationally expensive but accurate numerical routines often have to be replaced with fast and inaccurate methods, trading inference time for solution accuracy. This paper provides techniques to improve the quality of optimized control policies given a fixed computational budget. We achieve the above via a hypersolvers (Poli et al., 2020a) approach, which hybridizes a differential equation solver and a neural network. The performance is evaluated in direct and receding–horizon optimal control tasks, where the proposed approach shows consistent Pareto improvements in solution accuracy and control performance.

heading: 1 INTRODUCTION
text: Optimal control of complex, high–dimensional systems requires the implementation of computationally expensive numerical methods for differential equations (Pytlak, 2006; Rao, 2009). Here, real–time and hardware constraints preclude the use of accurate and expensive methods, forcing instead the application of cheaper and less accurate algorithms. While the paradigm of optimal control has successfully been applied in various domains (Vadali et al., 1999; Lewis et al., 2012; Zhang et al., 2016), improving accuracy while satisfying computational budget constraints is still a great challenge (Ross & Fahroo, 2006; Baotić et al., 2008).
To alleviate such computational limitations, we detail a procedure for offline optimization and subsequent online application of hypersolvers (Poli et al., 2020a) to optimal control problems. These hybrid solvers aim to achieve the accuracy of higher–order methods by augmenting numerical results of a base solver with a learning component trained to approximate local truncation residuals. When the overhead of a single forward–pass of the learning component is kept sufficiently small, hypersolvers improve the computation–accuracy Pareto front of low–order explicit solvers (Butcher, 1997). However, direct application of hybrid solvers to controlled dynamical system involves learning truncation residuals on the higher–dimensional space of state and control inputs.
To extend the range of applicability of hypersolvers to controlled dynamical systems, we propose two pretraining strategies designed to improve, in the set of admissible control inputs, on the average or worst–case hypersolver solution. With the proposed methodology, we empirically show that Pareto front improvements of hypersolvers hold even for optimal control tasks. In particular, we then carry out performance and generalization evaluations in direct and model predictive control tasks. Here, we confirm Pareto front improvements in terms of solution accuracy and subsequent control performance, leading to higher quality control policies and lower control losses.

heading: 2 NUMERICAL OPTIMAL CONTROL
text: We consider control of general nonlinear systems of the form
ẋ(t) = f(t, x(t), uθ(t)) x(0) = x0 (1)
with state x ∈ X ⊂ Rnx , input uθ ∈ U ⊂ Rnu defined on a compact time domain T := [t0, T ] where θ is a finite set of free parameters of the controller. Solutions of (1) are denoted with x(t) = Φ(x(s), s, t) for all s, t ∈ T . Given some objective function J : X × U → R; x0, uθ 7→ J(x0, uθ(t)) and a distribution ρ0(x0) of initial conditions with support in X , we consider the following nonlinear program, constrained to the system dynamics:
min uθ(t) Ex0∼ρ0(x0) [J(x0, uθ(t))]
subject to ẋ(t) = f(t, x(t), uθ(t)) x(0) = x0
t ∈ T
(2)
where the controller parameters θ are optimized. We will henceforth omit the subscript θ and write u(t) = uθ(t). Since analytic solutions of (2) exist only for limited classes of systems and objectives, numerical solvers are often applied to iteratively find a solution. For these reasons, problem 2 is often referred to as numerical optimal control.
Direct optimal control If the problem (2) is solved offline by directly optimizing over complete trajectories, we call it direct optimal control. The infinite–dimensional optimal control problem is time–discretized and solved numerically: the obtained control policy is then applied to the real target system without further optimization.
Model predictive control Also known in the literature as receding horizon control, Model Predictive Control (MPC) is a class of flexible control algorithms capable of taking into consideration constraints and nonlinearities (Mayne & Michalska, 1988; Garcia et al., 1989). MPC considers finite time windows which are then shifted forward in a receding manner. The control problem is then solved for each window by iteratively forward–propagating trajectories with numerical solvers i.e. predicting the set of future trajectories with a candidate controller u(t) and then adjusting it iteratively to optimize the cost function J (further details on the MPC formulation in Appendix B.2). The optimization is reiterated online until the end of the control time horizon.

heading: 2.1 SOLVER RESIDUALS
text: Given nominal solutions Φ of (1) we can define the residual of a numerical ODE solver as the normalized error accumulated in a single step size of the method, i.e.
Rk = R(tk, x(tk), u(tk)) = 1
p+1
[ Φ(x(tk), tk, tk+1)− x(tk)− ψ (tk, x(tk), u(tk)) ] (3)
where is the step size and p is the order of the numerical solver corresponding to ψ . From the definition of residual in (3), we can define the local truncation error ek := ∥∥ p+1Rk∥∥2 which is the error accumulated in a single step; while the global truncation error Ek = ‖x(tk)− xk‖2 represents the error accumulated in the first k steps of the numerical solution. Given a p–th order explicit solver, we have ek = O( p+1) and Ek = O( p) (Butcher, 1997).

heading: 3 HYPERSOLVERS FOR OPTIMAL CONTROL
text: We extend the range of applicability of hypersolvers (Poli et al., 2020a) to controlled dynamical systems. In this Section we discuss the proposed hypersolver architectures and pre–training strategies of the proposed hypersolver methodology for numerical optimal control of controlled dynamical systems.

heading: 3.1 HYPERSOLVERS
text: Given a p–order base solver update map ψ , the corresponding hypersolver is the discrete iteration
xk+1 = xk + ψ (tk, xk, uk)︸ ︷︷ ︸ base solver step + p+1 gω (tk, xk, uk)︸ ︷︷ ︸ approximator
(4)
where gω (tk, xk, uk) is some o(1) parametric function with free parameters ω. The core idea is to select gω as some function with universal approximation properties and fit the higher-order terms of the base solver by explicitly minimizing the residuals over a set of state and control input samples. This procedure leads to a reduction of the overall local truncation error ek, i.e. we can improve the base solver accuracy with the only computational overhead of evaluating the function gω . It is also proven that, if gω is a δ–approximator of R, i.e. ∀k ∈ N≤K
‖R (tk, x(tk), u(tk))− gω (tk, x(tk), u(tk))‖2 ≤ δ (5) then ek ≤ o(δ p+1), where δ > 0 depends on the hypersolver training results (Poli et al., 2020a, Theorem 1). This result practically guarantees that if gω is a good approximator for R, i.e. δ 1, then the overall local truncation error of the hypersolved ODE is significantly reduced with guaranteed upper bounds.

heading: 3.2 OPTIMAL CONTROL OPTIMIZATION WITH HYPERSOLVERS
text: Our approach relies on the pre–trained hypersolver model for obtaining solutions to the trajectories of the optimal control problem (2). After the initial training stage, control policies are numerically optimized to minimize the cost function J (see Appendix B.3 for further details). Figure 1 shows an overview of the proposed approach consisting in pre–training and system control.

heading: 4 HYPERSOLVER PRE–TRAINING AND ARCHITECTURES
text: We introduce in Section 4.1 loss functions which are used in the proposed pre–training methods of Section 4.2 and Section 4.3. We also check the generalization properties of hypersolvers with different architectures in Section 4.4. In Section 4.5 we introduce multi–stage hypersolvers in which an additional first–order learned term is employed for correcting errors in the vector field.

heading: 4.1 LOSS FUNCTIONS
text: Residual fitting Training the hypersolver on a single nominal trajectory {x(tk)}k results in a supervised learning problem where we minimize point–wise the Euclidean distance between the residual (3) and the output of gω , resulting in an optimization problem minimizing a loss function ` of the form
` (t, x, u) = 1
K K−1∑ k=0 ‖R (tk, x(tk), u(tk))− gω (tk, x(tk), u(tk))‖2 (6)
which is also called residual fitting since the target of gw is the residual R.
Trajectory fitting The optimization can also be carried out via trajectory fitting as following
` (t, x, u) = 1
K K−1∑ k=0 ‖x(tk+1), xk+1‖2 (7)
where x(tk+1) corresponds to the exact one–step trajectory and xk+1 is its approximation, derived via (4) for standard hypersolvers or via (12) for their multi–stage counterparts. This method can also be used to contain the global truncation error in the T domain. We will refer to ` as a loss function of either residual or trajectory fitting types; we note that these loss functions may also be combined depending on the application.
The goal is to train the hypersolver network to explore the state–control spaces so that it can effectively minimize the truncation error. We propose two methods with different purposes: stochastic exploration aiming at minimizing the average truncation error and active error minimization whose goal is to reduce the maximum error i.e., due to control inputs yielding high losses.

heading: 4.2 STOCHASTIC EXPLORATION
text: Stochastic exploration aims to minimize the average error of the visited state–controller space i.e., to solve the following optimization problem
min w E(x,u)∼ξ(x,u)` (t, x, u) (8)
where ξ(x, u) is a distribution with support in X × U of the state and controller spaces and ` is the training loss function. In order to guarantee sufficient exploration of the state–controller space, we use Monte Carlo sampling (Robert & Casella, 2013) from the given distribution. In particular, batches of initial conditions {xi0}, {ui0} are sampled from ξ and the loss function ` is calculated with the given system and step size . We then perform backpropagation for updating the parameters of the hypersolver using a stochastic gradient descent (SGD) algorithm e.g., Adam (Kingma & Ba, 2017) and repeat the procedure for every training epoch. Figure 2 shows pre–training results with stochastic exploration for different step sizes (see Appendix C.2). We notice how higher residual values correspond to higher absolute values of control inputs. Many systems in practice are subject to controls that are constrained in magnitude either due to physical limitations of the actuators or safety restraints of the workspace. This property allows us to design an exploration strategy that focuses on worst-case scenarios i.e. largest control inputs.

heading: 4.3 ACTIVE ERROR MINIMIZATION
text: The goal of active error minimization is to actively reduce the highest losses in terms of the control inputs, i.e., to solve the following optimization problem:
min ω max u∈U `(t, x, u) (9)
Similarly to stochastic exploration, we create distribution ξ(x, u) with support in X × U and perform Monte Carlo sampling of batches {xi0}, {ui0} from ξ. Then, losses are computed for each state xi0 with each control input u i 0 in order to obtain a new batch of controllers {ui′0 } yielding the maximum loss for each state, i.e. with elements:
ui ′
0 = arg max u0
`( , xi0, u0), u0 ∈ {ui0} (10)
This operation can be performed efficiently in parallel. After obtaining the new batch {ui′0 }, the loss is recalculated and SGD updates to hypersolver parameters are performed. Figure 3 shows a comparison of the pre–training techniques (further experimental details in Appendix C.2). The propagated error on trajectories for the hypersolver pre–trained via stochastic exploration is lower on average with random control inputs compared with the one pre– trained with active error minimization. The latter accumulates lower error for controllers yielding high residuals.

heading: 4.4 GENERALIZATION PROPERTIES OF DIFFERENT ARCHITECTURES
text: We have assumed the state and controller spaces to be bounded and that training be performed by sampling for their known distributions. While this is sufficient for optimal control problems given a priori known bounds, we also investigate how the system generalizes to unseen states and control input values. In particular, we found that activation functions have an impact on the end result of generalization beyond training boundaries.
We take into consideration two commonly used activation functions, Tanh : x→ ex−e−xex+e−x and ReLU : x → max(0, x), along with network architectures which employ activation functions containing periodic components: SIREN : x → sin (Wx+ b) (Sitzmann et al., 2020) and Snake : x → x + 1 a sin
2(ax) (Ziyin et al., 2020). We train hypersolver models with the different activation functions for the inverted pendulum model of (20) with common experimental settings (see Appendix C.3). Figure 4 and Figure 5 show generalization outside the training states, controllers and step sizes. We notice that while Tanh and ReLU perform well on the training set of interest, performance degrades rapidly outside of it. On the other one hand, SIREN and Snake manage to extrapolate the periodicity of the residual distribution even outside of the training region, thus providing further empirical evidence of the universal extrapolation theorem (Ziyin et al., 2020, Theorem 3).

heading: 4.5 MULTI–STAGE HYPERSOLVERS
text: We have so far considered the case in which the vector field (1) fully characterizes the system dynamics. However, if the model does not completely describe the actual system dynamics, first–order errors are introduced. We propose Multi-Stage Hypersolvers to correct these errors: an additional term is introduced in order to correct the inaccurate dynamics f . The resulting procedure is a modified version of (4) in which the base solver ψ (tk, xk, uk) does not iterate over the modeled vector field f but over its corrected version f?:
f? (tk, xk, uk) = f (tk, xk, uk)︸ ︷︷ ︸ Euler step +hw (tk, xk, uk)︸ ︷︷ ︸ first stage
(11)
where hw is a function with universal approximation properties. While the first stage hw is a first order error approximator, the second stage gω further reduces errors approximating the p–th order residuals:
xk+1 = xk + ψ tk, xk, uk, f? (tk, xk, uk)︸ ︷︷ ︸ corrected vector field + p+1 gω (tk, xk, uk)︸ ︷︷ ︸ second stage
(12)
We note that f? is continuously adjusted due to the optimization of hw. For this reason, it is not possible to derive the analytical expression of the residuals to train the stages with the residual fitting loss function (6). Instead, both stages can be optimized at the same time via backpropagation calculated on one–step trajectory fitting loss (7) which does not require explicit residuals calculation.

heading: 5 EXPERIMENTS
text: We introduce the experimental results divided for each system into hypersolver pre–training and subsequent optimal control. We use as accurate adaptive step–size solvers the Dormand/Prince method dopri5 (Dormand & Prince, 1980) and an improved version of it by Tsitouras tsit5 (Tsitouras, 2011) for training the hypersolvers and to test the control performance at runtime.

heading: 5.1 DIRECT OPTIMAL CONTROL OF A PENDULUM
text: Hypersolver pre–training We consider the inverted pendulum model with a torsional spring described in (20). We select ξ(x, u) as a uniform distribution with support in X × U where
X = [−2π, 2π] × [−2π, 2π] and U = [−5, 5] to guarantee sufficient exploration of the statecontroller space. Nominal solutions are calculated using tsit5 with absolute and relative tolerances set to 10−5. We train the hypersolver on local residuals via stochastic exploration using the Adam optimizer with learning rate of 3× 10−4 for 3× 105 epochs.
Direct optimal control The goal is to stabilize the inverted pendulum in the vertical position x? = [0, 0]. We choose t ∈ [0, 3] and a step size = 0.2 s for the experiment. The control input is assumed continuously time–varying. The neural controller is optimized via SGD with Adam with learning rate of 3 × 10−3 for 1000 epochs. Figure 6 shows nominal controlled trajectories of HyperEuler and other baseline fixed–step size solvers. Trajectories obtained with the controller optimized with HyperEuler reach final positions q = (1.6± 17.6)× 10−2 while Midpoint and RK4 ones q = (−0.6 ± 12.7) × 10−2 and q = (1.1 ± 12.8) × 10−2 respectively. On the other hand, the controller optimized with the Euler solver fails to control some trajectories obtaining a final q = (6.6± 19.4)× 10−1. HyperEuler considerably improved on the Euler baseline while requiring only 1.2% more Floating Point Operations (FLOPs) and 49.5% less compared to Midpoint. Further details are available in Appendix C.3.

heading: 5.2 MODEL PREDICTIVE CONTROL OF A CART-POLE SYSTEM
text: Hypersolver pre–training We consider an inaccurate version of the cart–pole system of (21) with a wrong parameter for the pole mass, hence we use the multi–stage hypersolver approach to correct the first–order error in the vector field. We select ξ(x, u) as a uniform distribution with support in X × U where X = [−2π, 2π] × [−2π, 2π] × [−2π, 2π] × [−2π, 2π] and U = [−10, 10]. Nominal solutions are calculated on the accurate system using tsit5 with absolute and relative tolerances set to 10−6. We train multi–stage HyperEuler (see Appendix A.2 for its explicit formulation) for the inaccurate system on nominal trajectories of the accurate system via stochastic exploration using the Adam optimizer with learning rate of 10−4 for 106 epochs.
Model predictive control The goal is to stabilize the cart–pole system in the vertical position around the origin, e.g. x? = [0, 0, 0, 0]. We choose t ∈ [0, 5] and a step size = 0.01 s for the experiment. The control input is assumed piece-wise constant during MPC sampling times. The receding horizon is chosen as 2 s. The neural controller is optimized via SGD with Adam with learning rate of 3 × 10−3 for 20 iterations at each sampling time. Figure 7 shows nominal controlled trajectories of multi–stage HyperEuler and other solvers. The adaptive–step tsit5 on the inaccurate system fails to stabilize the system at the origin position x = (−31.9 ± 8.0) cm, while multi–stage HyperEuler manages to stabilize the cart–pole system and improve on final positions x = (−5.3± 6.8) cm. Further details are available in Appendix C.4.

heading: 5.3 MODEL PREDICTIVE CONTROL OF A QUADCOPTER
text: Hypersolver pre–training We consider the quadcopter model of (23). We select ξ(x, u) as a uniform distribution with support in X × U where X is chosen as a distribution of possible visited states and each of the four motors i ∈ [0, 3] has control inputs ui ∈ [0, 2.17] × 105 rpm. Nominal solutions are calculated on the accurate system using dopri5 with relative and absolute tolerances set to 10−7 and 10−9 respectively. We train HyperEuler on local residuals via stochastic exploration using the Adam optimizer with learning rate of 10−3 for 105 epochs.
Model predictive control The control goal is to reach a final positions [x, y, z]? = [8, 8, 8] m. We choose t ∈ [0, 3] and a step size = 0.02 s for the experiment. The control input is assumed piece–wise constant during MPC sampling times. The receding horizon is chosen as 0.5 s. The neural controller is optimized via SGD with Adam with learning rate of 10−2 for 20 iterations at each sampling time. Figure 8 shows local residual distribution and control performance on the quadcopter over 30 experiments starting at random initial conditions which are kept common for the different ODE solvers. HyperEuler requires a single function evaluation per step as for the Euler solver compared to two function evaluations per step for Midpoint and four for RK4. Controlled trajectories optimized with Euler, Midpoint and RK4 collect an error on final positions of (1.09 ± 0.37) m, (0.71± 0.17) m, (0.70± 0.19) m respectively while HyperEuler achieves the lowest terminal error value of (0.66± 0.24) m. Additional experimental details are available in Appendix C.5.

heading: 5.4 BOUNDARY CONTROL OF A TIMOSHENKO BEAM
text: Hypersolver pre–training We consider the finite element discretization of the Timoshenko beam of (25). We create ξ(x, u) as a distribution with support in X × U which is generated at training time via random walks from known boundary conditions in order to guarantee both physical feasibility and sufficient exploration of the state-controller space (see Appendix C.6 for further details). Nominal solutions are calculated using tsit5 with absolute and relative tolerances set to 10−5. We train the hypersolver on local residuals via stochastic exploration using the Adam optimizer for 105 epochs, where we set the learning rate to 10−3 for the first 8× 104 epochs, then decrease it to 10−4 for 104 epochs and to 10−5 for the last 104.
Boundary direct optimal control The task is to stabilize the beam in the straight position, i.e. each of its elements i have velocities vit, v i r and displacements σ i t, σ i r equal to 0. We choose t ∈ [0, 3] and step size = 5 ms for the experiment. The control input is assumed continuously time– varying. The neural controller is optimized via SGD with Adam with learning rate of 10−3 for 1000 epochs. Figure 9 shows nominal controlled trajectories for HyperEuler and other baseline fixed– step size solvers. Control policies trained with Euler and Midpoint obtain averaged final states of (−2.8±4.2)×10−1 and (−0.04±4.6)×10−1 thus failing to stabilize the beam, while HyperEuler and RK4 obtain (−0.6±4.9)×10−3 and (−0.5±3.3)×10−3 respectively. HyperEuler considerably
improves on both the Euler and Midpoint baselines obtaining a very similar performance to RK4, while requiring 72.9% less FLOPs; the mean runtime per training iteration was cut from 8.24 s for RK4 to just 2.53 s for HyperEuler. Further details on this experiment are available in Appendix C.6.

heading: 6 RELATED WORK
text: This work is rooted in the broader literature on surrogate methods for speeding up simulations and solutions of dynamical systems (Grzeszczuk et al., 1998; James & Fatahalian, 2003; Gorissen et al., 2010). Differently from these approaches, we investigate a methodology to enable faster solution during a downstream, online optimization problem involving a potential mismatch compared to data seen during pre–training. We achieve this through the application of the hypersolver (Poli et al., 2020a) paradigm. Modeling mismatches between approximate and nominal models is explored in (Saveriano et al., 2017) where residual dynamics are learned efficiently along with the control policy while (Fisac et al., 2018; Taylor et al., 2019) model systems uncertainties in the context of safety–critical control. In contrast to previous work, we model uncertainties with the proposed multi–stage hypersolver approach by closely interacting with the underlying ODE base solvers and their residuals to improve solution accuracy.
The synergy between machine learning and optimal control continues a long line of research on introducing neural networks in optimal control (Hunt et al., 1992), applied to modeling (Lin & Cunningham, 1995), identification (Chu et al., 1990) or parametrization of the controller itself (Lin et al., 1991). Existing surrogate methods for systems (Grzeszczuk et al., 1998; James & Fatahalian, 2003) pay a computational cost upfront to accelerate downstream simulation. However, ensuring transfer from offline optimization to the online setting is still an open problem. In our approach, we investigate several strategies for an accurate offline–online transfer of a given hypersolver, depending on desiderata on its performance in terms of average residuals and error propagation on the online application. Beyond hypersolvers, our approach further leverages the latest advances in hardware and machine learning software (Paszke et al., 2019) by solving thousands of ODEs in parallel on graphics processing units (GPUs).

heading: 7 CONCLUSION
text: We presented a novel method for obtaining fast and accurate control policies. Hypersolver models were firstly pre–trained on distributions of states and controllers to approximate higher–order residuals of base fixed–step ODE solvers. The obtained models were then employed to improve the accuracy of trajectory solutions over which control policies were optimized. We verified that our method shows consistent improvements in the accuracy of ODE solutions and thus on the quality of control policies optimized through numerical solutions of the system. We envision the proposed approach to benefit the control field and robotics in both simulated and potentially real–world environments by efficiently solving high–dimensional space–continuous problems.

heading: CODE OF ETHICS
text: We acknowledge that all the authors of this work have read and commit to adhering to the ICLR Code of Ethics.

heading: REPRODUCIBILITY STATEMENT
text: We share the codes used in this paper publicly available as supplementary material for reproducibility. The following appendix also supplements the main text by providing additional clarifications. In particular, Appendix A provides further details on the considered hypersolver models. We provide additional information on optimal control policy in Appendix B while in Appendix C we provide details on on the system dynamics, architectures and other experimental details. Additional explanations are also provided in the submitted code.

heading: A ADDITIONAL HYPERSOLVER MATERIAL

heading: A.1 EXPLICIT HYPEREULER FORMULATION
text: Our analysis in the experiments takes into consideration the hypersolved version of the Euler scheme, namely HyperEuler. Since Euler is a first–order method, it requires the least number of function evaluations (NFE) of the vector field f in (1) and yields a second order local truncation error ek :=
∥∥ 2Rk∥∥2. This error is larger than other fixed–step solvers and thus has the most room for potential improvements. The base solver scheme ψ of (4) can be written as ψ (tk, xk, uk) = f (tk, xk, uk), which is approximating the next state by adding an evaluation of the vector field multiplied by the step size . We can write the HyperEuler update explicitly as
xk+1 = xk + f(tk, xk, uk) + 2gw (tk, xk, uk) (13)
while we write its residual as
R (x(tk), u(tk))) = 1
2 (Φ(x(tk), tk, tk+1)− x(tk)− f(tk, xk, uk)) (14)

heading: A.2 EXPLICIT MULTI–STAGE HYPEREULER FORMULATION
text: Similarly to the previous section, we show the explicit multi–stage HyperEuler formulation as
xk+1 = xk +  corrected vector field︷ ︸︸ ︷f (tk, xk, uk)︸ ︷︷ ︸ Euler step +hw (tk, xk, uk)︸ ︷︷ ︸ first stage + 2 gω (tk, xk, uk)︸ ︷︷ ︸ second stage
(15)
We also note that in this particular case the first stage hw and gω can be jointly trained; in particular, they can be chosen as outputs of a single network which is split in the two stages during the evaluation.

heading: A.3 HYPERSOLVERS FOR TIME–INVARIANT SYSTEMS
text: A time–invariant system with time–invariant controller can be described as following
ẋ(t) = f(x(t), u(x(t))) x(0) = x0 (16)
in which f and u do not explicitly depend on time. The models considered in the experiments satisfy this property.

heading: B CONTROL POLICY DETAILS

heading: B.1 OPTIMAL CONTROL COST FUNCTION
text: The general form of the integral cost functional can be written as follows
J(x(t), u(t)) = [x>(tf )− x?]P[x(tf )− x?] + ∫ tf t0 ( [x>(t)x?]Q[x(t)− x?] + u>(t)Ru(t) ) dt
(17) where matrix P is a penalty on deviations from the target x? of the last states, Q penalizes all deviations from the target of intermediate states while R is a regulator for the control inputs. Evaluation of (17) usually requires numerical solvers such as the proposed hypersolvers of this work. Discretizations of the cost functional are also called cost function in the literature.

heading: B.2 MODEL PREDICTIVE CONTROL FORMULATION
text: The following problem is solved online and iteratively until the end of the time span
min uk T−1∑ k=0 J (xk, uk)
subject to ẋ(t) = f(t, x(t), u(t)) x(0) = x0
t ∈ T
(18)
where J is a cost function and T ∈ T is the receding horizon over which predicted future trajectories are optimized.

heading: B.3 NEURAL OPTIMAL CONTROL
text: We parametrize the control policy of problem (2) as uθ : t, x 7→ uθ(t, x) where θ is a finite set of free parameters. Specifically, we consider the case of neural optimal control in which controller uθ is a multi–layer perceptron. The optimal control task is to minimize the cost function J described in (17) and we do so by optimizing the parameters θ via SGD; in particular, we use the Adam (Kingma & Ba, 2017) optimizer for all the experiments.

heading: C EXPERIMENTAL DETAILS
text: In this section we include additional modeling and experimental details divided into the different dynamical systems.

heading: C.1 HYPERSOLVER NETWORK ARCHITECTURE
text: We design the hypersolver networks gw as feed–forward neural networks. Table 1 summarizes the parameters used for the considered controlled systems, where Activation denotes the activation functions, i.e. SoftPlus: x 7→ log(1 + ex), Tanh: x 7→ ex−e−xex+e−x and Snake : x→ x+ 1a sin2(ax) (Ziyin et al., 2020).
We also used the vector field f as an input of the hypersolver, which does not require a further evaluation since it is pre–evaluated at runtime by the base solver ψ. We emphasize that the size of the network should depend on the application: a too–large neural network may require more computations than just increasing the numerical solver’s order: Pareto optimality of hypersolvers also depends on their complexity. Keeping their neural network small enough guarantees that evaluating the hypersolvers is cheaper than resorting to more complex numerical routines.
1This architecture refers to the optimal control experiment. Details on hypersolver models for the generalization experiment on the inverted pendulum are available in Appendix C.3.
2The output layer considers both the first stage hw and the second stage gω of the multi–stage HyperEuler scheme which can be trained as a single neural network whose outputs are split while solving the ODE problem.

heading: C.2 SPRING-MASS SYSTEM
text: System Dynamics The spring-mass system considered is described in the Hamiltonian formulation by [
q̇ ṗ
] = [ 0 1/m −k 0 ] [ q p ] + [ 0 1 ] u (19)
where m = 1 [Kg] and k = 0.5 [N/m].
Pre-training methods comparison We select ξ(x, u) as a uniform distribution with support in X ×U where X = [−20, 20]× [−20, 20] while U = [−100, 100]. Nominal solutions are calculated on the accurate system using dopri5 with relative and absolute tolerances set to 10−7 and 10−9 respectively. We train two separate HyperEuler models with different training methods on local residual for step size = 0.03 s: stochastic exploration and active error minimization. The optimizer used is Adam with learning rate of 10−3 for 104 epochs.
Hypersolvers with different step sizes We select ξ(x, u) as a uniform distribution with support in X × U where X = [−5, 5] × [−5, 5] while U = [−20, 20]. Nominal solutions are calculated on the accurate system using dopri5 with relative and absolute tolerances set to 10−7 and 10−9 respectively. We train separate HyperEuler models with stochastic exploration with different step sizes . The optimizer used is Adam with learning rate of 10−3 for 104 epochs.
C.3 INVERTED PENDULUM
System Dynamics We model the inverted pendulum with elastic joint with Hamiltonian dynamics via the following: [
q̇ ṗ
] = [ 0 1/m −k −β/m ] [ q p ] − [ 0 mgl sin q ] + [ 0 1 ] u (20)
where m = 1 [Kg], k = 0.5 [N/ rad], r = 1 [m], β = 0.01 [Ns/ rad], g = 9.81 [m/s2].
Pre–training for the generalization study We perform sampling via stochastic exploration from the uniform distribution ξ(x, u) with support in X × U with X = [−2π, 2π] × [−2π, 2π] and U = [−10, 10] for the different architectures. We choose as a common time step = 0.1 s; the networks are trained for 100000 epochs with the Adam optimizer and learning rate of 10−3. The network architectures share the same parameters as the inverted pendulum ones in 1, while the activation functions are substituted by the ones in Figure 4. The SIREN architecture is chosen with 2 hidden layers of size 64.
Additional visualization Figure 10 provides an additional visualization of the inverted pendulum controlled trajectories from Figure 6 with positions q and momenta p over time.

heading: C.4 CART-POLE
text: System Dynamics We consider a continuous version of the cart-pole system in OpenAI Gym (Brockman et al., 2016), appropriately modified for PyTorch and batched training. The dynamic
model is described by the four variables x, ẋ, θ, θ̇ (Florian, 2005) and the accelerations update is as following
θ̈ = g sin θ + cos θ
( −F−mplθ̇2 sin θ
mc+mp ) l (
4 3 − mp cos2 θ mc+mp ) ẍ = F +mpl ( θ̇2 sin θ − θ̈ cos θ ) mc +mp
(21)
where mc = 1 [Kg], mp = 0.1 [Kg], l = 0.5 [m] and g = 9.81 [m/s2].
Additional experimental details For the multi–stage HyperEuler experiment, we increase the mass of the pole mp to 1 [Kg] for the nominal system while the hypersolver only has knowledge of the mass mp = 0.1 [Kg]. We also pre-train a single neural network whose first 4 outputs constitute the first stage hw and the last 4 are gω . Figure 11 shows mean local residuals for the learned first stage only and for the complete multi–stage scheme.
Comparison with other baselines We additionally compare the multi–stage hypersolver approach on the inaccurate model with other baseline methods, namely the adaptive tsit5with knowledge of the inaccurate model only, the Euler method iterating on learned system residual dynamics, and the Euler and Midpoint solvers with the accurate (nominal) model. For the residual dynamics approach, the goal is to learn the corrected vector field f?. We parametrize the residual dynamics as a neural network with the same number of layers, hidden size and same activation functions as multi–stage HyperEuler described in Table 1. We optimize the network via SGD with the Adam optimizer until convergence by approximating the vector field directly by drawing batches from the same distribution ξ(x, u) on which multi–stage HyperEuler is trained, i.e. by minimizing the loss function
`(t, x, u) = 1
K K−1∑ k=0 ‖ψ (tk, xk, uk, f?(tk, xk, uk)) , ψ (tk, xk, uk, fnom(tk, xk, uk))‖2 (22)
in which the corrected vector field f?(·) = f(·) + rw(·) with f being the inaccurate dynamics and rw the residual dynamics network we optimize, fnom is the nominal dynamics and ψ is the solution corresponding to the base solver tsit5 set with absolute and relative tolerances of 10−5. Figure 12 shows error propagation over time of the different approaches. While propagating the corrected vector field learning the residual dynamics and multi–stage HyperEuler both improve on the baseline tsit5 solver iterating on the inaccurate model, we note that multi–stage HyperEuler manages to also outperform the Euler solver on the accurate dynamics.

heading: C.5 QUADCOPTER
text: System Dynamics The quadcopter model is a suitably modified version of the explicit dynamics update in (Panerati et al., 2021) for batched training in PyTorch. The following accelerations update describes the dynamic model
ẍ = ( R · [0, 0, kF ∑3 i=0ω 2 i ]− [0, 0,mg] ) m−1
ψ̈ = J−1 ( τ(l, kF , kT , [ω 2 0 , ω 2 1 , ω 2 2 , ω 2 3 ])− ψ̇ × ( Jψ̇ )) (23)
mwhere x = [x, y, z] corresponds to the drone positions and ψ = [φ, θ, ψ] to its angular positions; R and J are its rotation and inertial matrices respectively, τ(·) is a function calculating the torques induced by the motor speeds ωi, while arm length l, mass m, gravity acceleration constant g along with kF and kT are scalar variables describing the quadcopter’s physical properties.

heading: C.6 TIMOSHENKO BEAM
text: System Dynamics We consider as a system from the theory of continuum dynamics the Timoshenko beam with no dissipation described in (Macchelli & Melchiorri, 2004). The system can be described in the coenergy formulation by the following partial differential equation (PDE)ρA 0 0 00 Iρ 0 00 0 Cb 0
0 0 0 Cs
 ∂ ∂t vtvrσr σt  =  0 0 0 ∂x0 0 ∂x 10 ∂x 0 0 ∂x −1 0 0  vtvrσr σt  (24) where ρ is the mass density, A is the cross section area, Iρ is the rotational inertia, Cs and Cb are the shear and bending compliance; the discretized state variables vt, vr represent the translational and rotational velocities respectively while σt, σr denote the translational and rotational displacements.
In order to discretize the problem, we implement a software routine based on the fenics (Alnæs et al., 2015) open–source software suite to obtain the finite–elements discretization of the Timoshenko PDE of (24) given the number of elements, physical parameters of the model and initial conditions of the beam. We choose a 40 elements discretization of the PDE for a total of 160 dimensions of the discretized state z = [vt, vr, σt, σr]> and we initialize the beam at time t = 0 as z(x, 0) = [sin(πx), sin(3πx), 0, 0]. The system can thus be reduced to the following controlled linear system
 v̇tv̇rσ̇t σ̇r  =  × × × −M−1ρAD>1 × × −M−1Iρ D>2 −M −1 Iρ D>0
× M−1Cb D2 × × M−1Cs D1 M −1 Cs D0 × ×
 vtvrσt σr +  × M−1ρABF M−1Iρ BT × × × × × [u1∂u2∂ ]
(25) where the mass matrices MρA, MIρ , MCb , MCs , matrices D0, D1, D2, vectors BF , BT are computed through the fenics routine and boundary controllers u1∂ and u 2 ∂ are the control torque and the control force applied at the free end of the beam.
Stochastic exploration strategy via random walks We pre–train the hypersolver model via stochastic exploration of the state-controller space X × U . We restrict the boundary control input values u = [u1∂ , u 2 ∂ ] > in [−1, 1] × [−1, 1]. As for the state space, naively generating a probability distribution with box boundaries on each of the 160 dimensions of X would require an inefficient search over this high-dimensional space: in fact, not every combination is physically feasible due to the Timoshenko beam’s structure. We solve this problem by propagating batched trajectories with RK4 from the initial boundary condition z(x, 0) = [sin(πx), sin(3πx), 0, 0] with random control actions sampled from a uniform distribution with support in [−1, 1] × [−1, 1] applied for a time t1 ∼ U [0.002, 1] s. We save the states {z(x, t1)i} and forward propagate from these states again by sampling from the controller and time distributions. We repeat the process K times and obtain a sequence [{z(x, t1)i}, . . . , {z(x, tK)i}] of batched initial conditions characterized by physical feasibility. Finally, we train the hypersolver with stochastic exploration by sampling from the generated distribution ξ(x, u) on local one–step residuals as described in Section 5.4. This initial state generation strategy is repeated every 100 epochs for guaranteeing an extensive exploration of all possible boundary conditions. Figure 13 shows the error propagation over controlled trajectories: the trained HyperEuler achieves the lowest error among baseline fixed–step solvers.
Additional details on the results We report additional details regarding the runtime of the experiments on the Timoshenko Beam. As for the training time, training the hypersolver for 105 epochs takes around 80 minutes, where the time for each training epoch slightly varies depending on the length of the sequence of initial condition batches obtained via random walks. As for the averaged runtime per training epoch during the control policy optimization, HyperEuler takes (2.53± 0.09) s per training iteration, Euler (2.01 ± 0.04) s, Midpoint (4.02 ± 0.08) s and RK4 (8.24 ± 0.14) s. Experiments were run on the CPU of the machine described in Section C.7.

heading: C.7 HARDWARE AND SOFTWARE
text: Experiments were carried out on a machine equipped with an AMD RYZEN THREADRIPPER 3960X CPU with 48 threads and two NVIDIA RTX 3090 graphic cards. Software–wise, we used PyTorch (Paszke et al., 2019) for deep learning and the torchdyn (Poli et al., 2020b) and torchdiffeq (Chen et al., 2019) libraries for ODE solvers.

Review summaries:
In this paper, the authors apply the idea of hypersolvers to numerical optimal control. The idea is to approximate the dynamics with a low-order ODE solver and learn the truncation residual dynamics using a neural net. They also propose a multi-stage hypersolvers method to both learn the system residual dynamics (due to erroneous system model) and truncation residual dynamics (due to ODE solver). The performance of the proposed approach is shown in simulation experiments with a pendulum, a cart-pole, and a quadcopter with comparisons to ODE solvers of different orders.

This paper extends the hypersolvers framework (Poli et al., 2020a) to the setting of optimal control. The hypersolvers framework combines a cheap numerical dynamical system solver with a neural network trained to approximate the truncation error of this solver, which together yield a cheap and relatively accurate solver. This paper shows how to utilize the hypersolvers framework within a controlled dynamical system by exploring different choices of loss functions, training methods, and network architectures for hypersolver pre-training. The authors further provide an augmentation to this framework, called multi-stage hypersolvers, to further account for misspecifications of the dynamics. Via experiments, the authors then demonstrate the performance of optimal control strategies (direct optimal control and MPC) utilizing hypersolvers in different settings (notably the pendulum, cart-pole, and quadcopter settings). They show e.g. that even in cases where their method receives misspecified dynamics, they are able to achieve control performance that is comparable to that obtained in cases where the underlying solver had knowledge of the correct dynamics.

The paper descibes the application of hypersolves, which warm-start ODE solvers with neural nets that have distilled many precomputed solutions, to optimal control ODEs.  This warm-starting improves downstream accuracy when the computational budget is fixed.  This method is applied for trajectory optimziation for pendulum swing-up and MPC for cartpole and quadcopter systems.