Neural Solvers for Fast and Accurate Numerical Optimal Control
Authors: Federico Berto
Stefano Massaroli
Michael Poli
Jinkyoo Park
Venue: ICLR 2022 Poster (2022)
Keywords: Deep Learning
Numerical Methods
Optimal Control

Abstract: Synthesizing optimal controllers for dynamical systems often involves solving optimization problems with hard real-time constraints. These constraints determine the class of numerical methods that can be applied: computationally expensive but accurate numerical routines are replaced by 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 approach, which hybridizes a differential equation solver and a neural network. The performance is evaluated in direct and receding-horizon optimal control tasks in both low and high dimensions, where the proposed approach shows consistent Pareto improvements in solution accuracy and control performance.

Published as a conference paper at ICLR 2022
NEURAL SOLVERS FOR FAST AND ACCURATE NUMER- ICAL OPTIMAL CONTROL
Federico Berto
Stefano Massaroli
KAIST, Di fifffMML fberto0kai wtw.eokkkistt....... at.
The University of Tokyo, Di ffFqML. massaro1 eesaroll@roboer Br.bot.t. w.t.bot...ttokko.ac.....
Michael Poli
Jinkyoo Park KAIST
Stanford University, Di biffEo fEMML zymrae18cs stt.eeleesssss.forooooooooo edu
ABSTRACT
Synthesizing optimal controllers for dynamical systems often involves solving op- timization problems with hard real-time constraints. These constraints determine the class of numerical methods that can be applied: computationally expensive but accurate numerical routines are replaced by fast and inaccurate methods, trad- ing inference time for solution accuracy. This paper provides techniques to im- prove the quality of optimized control policies given fixed computational bud- get. We achieve the above via a hypersolvers (Poli et al., 2020a) approach, which hybridizes differential equation solver and neural network. The performance is evaluated in direct and saceding-horizz optimal control tasks in both low and high dimensions, where the proposed approach shows consistent Pareto improve- ments in solution accuracy and control performance.
INTRODUCTION
Optimal control of complex, high-dimenniooal systems requires computationally expensive numer- ical 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 great challenge (Ross & Fahroo, 2006; Baotić et al., 2008). To alleviate computational overheads, we detail a procedure for offline optimization and subsequent online application of hypersolvers (Poli et al., 2020a) to optimal con- trol problems. These hybrid solvers achieve the accuracy of higher-orde methods by augmenting numerical results of base solver with learning component trained approximate local truncation residuals. When the cost of a single forward-pass of the learning component is kept sufficiently small, hypersolvers improve the computation-aacuracc Pareto front of low-order explicit solvers (Butcher, 1997). However, direct application of hybrid solvers to controlled dynamical system in- volves learning truncation residuals on the higher- dimensionll spaces 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 per- formance, leading to higher quality control policies and lower control losses. In dineeninmsii regimes, we obtain the same control policy as the one obtained by accurate high-ordes solvers with more than 3x speedup.

---

Published as a conference paper at ICLR 2022
Hypersolver Training
Optimal Control
Dastibution (Gx.) of possible stares and controller values Rygensdvis Su waa base scvvie
Jix, s)
Initaal condition Xa
Sumple batches from giz,a]
Solve optiniinacion ion probllm
Pard-forward hypenslorr
I Mrs sal (0,7-1)
Tre tataco
Fast %
Ug
Calenlate E(R (x.))) and Lackgrooaggte
Apply contrul impure
Shift recaaling Furward
Calerlae residnale
Q
-
Figure I: Overview of the proposed method. [Left] The hypersolver is trained to approximate residuals given distribution of control inputs and states. [Right] The pre trained hypersolver model then used to accelerate improve the accuracy numerical solutions used during optimization of control policies, leading to higher- quality controllers.
2 NUMERICAL OPTIMAL CONTROL
We consider control of general nonlinear systems of the form
i(1) (t,), ua()))
x (0) 30
with state x E X C R"x input Ug E U C TT defined compact time domain T = [4o 7] where 0 is a finite set of free parameters of the controller. Solutions of (1) are denoted with xt)) = (<x(s), s, t) for all s. T. Given some objective function J X xul - R; 2o. ug ( ug (t)) and distribution Po (2o) of initial conditions with support in X, we consider the following nonlinear program, constrained to the system dynamics:
min (t)
E E [x un (())
subject to
=(t) = ftt xt), ut())
x00)
E T
where the controller parameters 0 are optimized. We will henceforth omit the subscript 0 and write u(t) (t).). Since analytic solutions of (2) exist only for limited classes of systems and objectives, numerical solvers often applied to iteratively find solution. For these reasons, problem is often referred to as numerical optimal control.
Direct optimal control the problem (2) is solved offline by directly optimizing over complete trajectories, we call it direct optimal control. The infinite-din fimniteddimessionll optimal control problem is disc-discet 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 Pre- dictive Control (MPC) is class of flexible control algorithms capable of taking into consideration constraints and nonlinearities (Mayne Michalska, 1988; Garcia et al., 1989). MPC considers fi- nite time windows which are then shifted forward in receding manner. The control problem is then solved for each window by iteratively forward -wwaarddrrppgatinng trajectories with numerical solvers iec. predicting the set of future trajectories with candidate controller u(t) and then adjusting it it- eratively 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.
2

---

Published as a conference paper at ICLR 2022
2.1 SOLVER RESIDUALS
Given nominal solutions & of (1) we can define the residual of numerical ODE solver as the normalized error accumulated in single step size of the method, i.e.
_ (<<(t+)) las (x+1) x(4x) (1a, s(()), (((k),
RR R((a, =044); u((+))
(3)
c-T
where e is the step size and is the order of the numerical solver corresponding to the From the which the single step; while the global truncation error [x((x) I& ll represents
definition of residual in (3), we can define the local truncation envor ek := c error accumulated the enror accumulated in the first k steps of the numerical solution. Given a p-th order explicit solver, we have eg O(c*+1) and & O(c)) (Butcher, 1997).
HYPERSOLVERS FOR OPTIMAL CONTROL
We extend the range of applicabilitt of hypersolvers (Poli et al., 2020a) controlled dynamical sys- tems. 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.
3.1 HYPERSOLVERS
Given p-order base solver update map t/., the corresponding hypersolver is the discrete iteration
Fk++1 = Ft. + ea/e (4R, Ih: us) +e Ii (J,, uk)
base solver step
approximator
where Aw DK, Tk, ux) is (11) parametric function with free parameters The core idea is select Ou as some function with universal approximation properties and fit the higher-order terms of the base solver by explicitly minimizing the residuals over set of state and control input samples. This procedure leads to a reduction of the overall local truncation error eg, i.c. we can improve the base solver accuracy with the only computational overhead of evaluating the function Jw- It is also proven that, if 9w is a 8 apapoxximator of R. i.e. Vk E N
|| (tx, r(t2), u(tx)) gw (tk, 2(tx), u(tx))ll2 <
(5)
then S 1() where > O depends on the hypersolver training results (Poli et al., 2020a, Theorem 1). This result practically guarantees that if gw is a good approximator for i.e. 8 < 1, then the overall local truncation error of the hypersolved ODE is significantly reduced with guaran- teed upper bounds.
3.2 NUMERICAL OPTIMAL CONTROL WITH HYPERSOLVERS
Our approach relies 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 I shows an overview the proposed approach consisting
training and system control.
HYPERSOLVER PRE-TRAINING AND ARCHITECTURES
We introduce in Section 4.1 loss functions which are used in the proposed pre ge-ttrining methods of Section 4.2 and Section 4.3. We also check the generalization properties of hypersolvers with different architectures Section 4.4. In Section 4.5 we introduce multi-stage hypersolvers which an additional first-order learned term is employed for correcting errors in the vector field.
4.1 LOSS FUNCTIONS
Residual fitting Training the hypersolver on a single nominal trajectory (=(tt)}e results a supervised learning problem where we minimize point-wise the Fuclidean distance between the
3

---

Published as a conference paper at ICLR 2022
0.01 [8)
0.03 [s]
0.1
0.3
1.0 [8]
10
Moan Rosichak X
Hyperbuier
10-"-
Ealer
Midpoint
RK4
[M]
[N]
Figure 2: Mcan local residuals of the spring-mass system of (17) as a function of control inputs at different step sizes HyperFuler (see Appendix A.1 for its explicit formulation) improves the local residuals compared the baseline Fuler and compared to higher-order ODE solvers larger step sizes.
residual (3) and the output of Iw. resulting in optimization problem minimizing loss function & of the form
1 K
e (t, u)
Σ (ta, (t2,, x(tx), u(()) Iw (t&, (txx), (tt)) ll2)))
which is also called residual fitting since the target of Iw is the residual R.
Trajectory fitting
The optimization can also be carried out via trajectory fitting as following
K
e (f, xr
1
Σ
(z=(t Z.
K
where x((x+1) corresponds to the exact one-step trajectory and 2x+1 is its approximation, derived via (4) for standard hypersolvers or via (11) 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 loss function of either residual or trajectory fitting types; we note that these loss functions may also be combined depending the application. The goal to train the hypersolver network explore the state- control spaces so that it can effectively minimize the truncation eπror. 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.
4.2 STOCHASTIC EXPLORATION
Stochastic exploration aims to minimize the average error of the visited state controlene space i.e., to produce optimal hypersolver parameters w" as the solution of nonlinear program
arg min Eg [e (t, x, u)]
where {xx, u) is 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 {2(), {ur} are sampled from & and the loss function E is calculated with the given system and step size e. We then perform backpropagation for updating the parameters of the hypersolver using stochastic gradien! descent (SGD) algorithm e.g., Adam (Kingma & Ba, 2017) and repeat the procedure for every training epoch. Figure 2 shows pre- pre-training results with stochastic cxploration for different step sizes (see Appendix C.2). We notice how higher residual values generally correspond to higher absolute values of control inputs. Many systems in practice subject to controls that constrained in magnitude either due 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.
4

---

Published as a conference paper at ICLR 2022
Without Hypersolver
Tanh
ReLU max(0,x)
SIREN (Wxx
Snake (ax)
Ealer Residwals
75
10
Mean Residuni
60
45
-10
15
4
2rt 4rt
1x
2x
1x
2
Figure 4: Generalization outside of the training region (red rectangle) in the state space of inverted pendu- lum model with different hypersolver activation functions. Architectures containing activation functions with periodic components achieve better extrapolation properties compared to others.
4.3 ACTIVE ERROR MINIMIZATION
The goal of active error minimization is to ac- tively reduce the highest losses in terms of the control inputs, i.e., to obtain w as the solution to minmax problem:
7.5
MAE
w0 arg min E(t, x.,
5.0
(9)
well
2.5
Similarly to stochastic exploration, we create distribution E(x, u) with support X' ll and perform Monte Carlo sampling of n batches {<=',""), from ç. Then, losses are computed pair-wisely for each state x! j O,..., with each control input u²,
0.0
x10
7.5
MAE
5.0
0......1. We then take the first controllers (ug) yield- ing the maximum loss for each state. The loss is recalculated using these controller values with their respective states and SGD updates to hypersolver parameters are performed. Fig- shows comparison of the pre- training techniques (further experimental details in Ap- pendix C.2). The propagated error on trajecto- ries for the hypersolver pre-traaiee treited via stochas- tic exploration is lower on average with ran- dom control inputs compared with the one pre- trained with active error minimization. The latter accumulates lower error for controllers yielding high residuals.
2.5
0.0
0.0
0.3
Time [s]
Stochastic Exploration
Active Ervor Minimization
Midpoint
Mean Absolute Error (MAE) along tra- jectories with different pre- trainiin techniques the spring-mas, mans system (17). (TTop] Stochastic explo- ration performs better on average i.c. u e [-100, 100). [Bottom] Active error minimization achieves better re- sults in limit situations in the case of bang-bang controller i.e. (-100, 100), in which controllers yielding highest residuals have been minimized.
Figure 3:
Different exploration strategies may he used depending on the down-stream control task.
4.4 GENERALIZATION PROPERTIES OF DIFFERENT ARCHITECTURES
We have assumed the state and controller spaces 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 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 I - c"-e-- architectures which employ activation functions containing periodic components: SIREN x - sin sin? (ax) (Ziyin et al., 2020). We train hypersolver models with the different activation functions for the inverted pendulum model of (18) with common experimental settings (see Appendix C.3). Figure 4 shows generalization
c=+e-= and ReLU x - max(0, 2),, along with network
sin (Waz b) (Sitzmann al., 2020) and Snake
5

---

Published as a conference paper at ICLR 2022
outside the training states (see Figure 9 in Appendix C.3 for generalization of 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 al., 2020, Theorem, 3).
Activation function choice plays an important role Hypersolver performance and generalization.
4.5 MULTI-STAGE HYPERSOLVERS
We have so far considered the case in which the vector field (1) fully characterizes the system dy- namics. However, the model does not completely describe the actual system dynamics, first-order emrors are introduced. We propose Multi- SelltiiSe Stagg Hypersolvers additional term is introduced order to correct the inaccurate dynamics f. The resulting procedure is modified version of (4) in which the base solver t/e (tx, Tk, uy) does not iterate over the modeled vector field f but over its corrected version f":
correct these errors: an
I (Lk, Tkk ux) = (1s, us.) + hio (1x, W.
(10)
partial dynamics
inner stage
where Aw is function with universal approximation properties. While the inner stage hw is first- order error approximator, the outer stage gw further reduces errors approximating the p-th order residuals:
Fk+1 Uk + ea/e (x: Tky Uk
S* [0 (F uk)
e
Iu (T uge
(11)
corrected dynamics
stage
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 be optimized at the same time via backpropagation calcu- lated on one- trajectory fitting loss (7) which does not require explicit residuals calculation.
5 EXPERIMENTS
We introduce the experimental results divided for each system into hypersolver pre-training and sub- sequent optimal control. We 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.
5.1 DIRECT OPTIMAL CONTROL OF PENDULUM
Hypersolver pre-training We consider the inverted pendulum model with torsional spring described in (18). We select E(x, u) as a uniform distribution with support in U where X' 1-2x, 2x] [-2x,2 2nr and U [-5,5] to guarantee sufficient exploration of the state- controller space. Nominal solutions are calculated using tsit5 with absolute and relative tolerances set to 10 We train the hypersolver
local residuals via stochastic exploration using the Adam optimizer with learning rate of 3 x 10-4 for 3 x 10a epochs.
Direct optimal control The goal is to stabilize the inverted pendulum in the vertical position x² [o. 0). We choose E [0.3] and step size 0.2 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 5 shows nominal controlled trajectories of HyperEuler and other baseline fixed-step step size solvers. Trajectories obtained with the controller optimized with HyperEuler reach final positions 9 (1.6 + 17.6) x 10 while Midpoint and RK4 ones 9 (-0.6 = 12.7) x 10 10-2 and = (1.1 + 12.8) x 10 respectively. On the other hand, the controller optimized with the Fuler solver fails control some trajectories obtaining final

---

Published as a conference paper at ICLR 2022
q q= (6.6 10-1 HyperEuler considerably improved the Euler baseline while requiring only 1.2% more Floating Point Operations (FLOPs) and 49.5% less compared to Midpoint. Further details available in Appendix C.3.
HyperEusler
Euler
Midpoint
RK4
Figure 5: Direct optimal control of the inverted pendulum in phase space. While the controller optimized with the Fuler solver fails to control the system for some trajectories, the obtained with HyperFuler improves the performance while introducing minimal overhead with results comparable higher-order solvers.
5.2 MODEL PREDICTIVE CONTROL OF A CART-POLE SYSTEM
Hypersolver pre-training We consider the partial dynamics of the cart-pole system of (19) with wrong parameters for the frictions between cart and track as well as the one between cart and pole. We employ the multi-stage hypersolver approach to correct the first-order eπrоr in the vector field as well as base solver residual. We select (ax, u) as uniform distribution with support in x 2 where X' = [-2 2n)] x 2m, 2x) x 1T, 2x 2m, 2T] and 4 = [-10. 10). Nominal solutions are calculated on the accurate system using RungeKutta instead of adaptive-step solvers due faster training times. We train our multi-stage Hypersolver (i.e. multi-stage hypersolver with second-order Midpoint as base solver with the partial dynamics) on nominal trajectories of accurate system via stochastic exploration using the Adam optimizer for 5 10 epochs, where set the learning rate to 10 for the first 10t epochs, then decrease it to for 10 epochs and 10-4 for the last 10
Model predictive control The goal is stabilize the cart-pole system in the vertical position around the origin, e.g. x² [0. 0, O, 0]. We choose t E [0,3] and step size 0.05 for the experiment. The control input assumed piece-wise constant during MPC sampling times. The receding horizon is chosen as 1 The neural controller is optimized via SGD with Adam with learning rate of 3 x 10-3 3 for maximum of 200 iterations at each sampling time. Figure shows nominal controlled trajectories of multi-stage Hypersolver and other baseline solvers. The Midpoint solver on the inaccurate model fails to stabilize the system at the origin position x (39.7 + 97.7) em, while multi-stage Hypersolver manages to stabilize the cart-pole system and improve on final positions x = (7.8 + 3.0) cm. Further details are available in Appendix C.4.
Imput [N]
1o0 mang mmnen
100
Time [c]
Time [g]
Time Ed
1iose [s]
Hypenselver (Inaccumbee
Midpore mudel)
Figure 6: Model Predictive Control with constrained inputs on the cart-pole model. MPC with the Midpoint solver iterating the partial dynamic model successfully swings up the pole but fails reach the target position. Multi-stage Hypersollvrr with the Midpoint base solver has knowledge restricted to the inaccurate system, yet it manages to obtain similar control performance compared to controllers with access nominal dynamics while also needing less control effort and absolute energy inflow compared to its base solver.
Multi-stage Hypersolvers can correct first-order error on dynamic models base solver residuals.

---

Published as a conference paper at ICLR 2022
5.3 MODEL PREDICTIVE CONTROL A QUADCOPTER
Hypersolver pre-training We consider the quadcopter model of (20), We select E(x, u) as a uniform distribution with support in X' U where X' chosen distribution of possible visited states and each of the four motors i E [0, 3] has control inputs u² E [0, 2.17 x 10 rpm. Nominal solutions are calculated on the accurate system using dopri5 with relative and absolute tolerances set to 10- and 10-9 respectively. We train Hyperfuler on local residuals via stochastic exploration using the Adam optimizer with learning rate of 10-3 for 10 epochs.
Model predictive control The control goal is to reach final positions [xr, y, 2|* [8,8, 8, We choose 1. E [0,3] and a step size = 0.02 for the experiment. The control input is assumed piece-wise constant during MPC sampling times. The receding horizon is chosen 0.5 The neural controller is optimized via SGD with Adam with learning rate of 10 2 for 20 iterations at each sampling time. Figure 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 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) (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.
2x10
Nain K
4x
S3331
HoywWBske
Galer
Milpors
RKC
Poyrefidkrr
Milpeie
3ki
Figure 7: [Left] Local residual distribution for the quadcopter model for 0.02 [Center] Trajectories controlled quadcopters with MPC whose receding horizon controller is optimized by solving the ODE with different methods. [Right] Final positions error distribution. The proposed approach with HyperEuler achieves lower average error compared to other baseline solvers while requiring a low overhead compared to higher- order solvers due to smaller number of dynamics function evaluations.
5.4 BOUNDARY CONTROL OF A TIMOSHENKO BEAM
Hypersolver pre-training We consider the finite element discretization of the Timoshenko beam of (22). We create E(x, u) as distribution with support in X' x U which is generated at training time via random walks from known boundary conditions in order guarantee both physical feasi- bility 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 10 epochs, where we set the learning rate to 10 3 for the first 8 x 1 epochs, then decrease it to 10 for 10a epochs and to 10 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 have velocities v.. of and displacements of of equal to O. We choose E [0.3] and step size e 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 8 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) x 10 and 0.004.44440 4.6) x 10-1 thus failing to stabilize the beam, while HyperEuler and RK4 obtain (-0.6=4.9) x 10 and (-0.5+3.3) x 10 respectively. HyperEuler considerably improves on both the Euler and Midpoint baselines obtaining 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 just 2.53 for Hyperfuler. Further details this experiment available in Appendix C.6.
8

---

Published as a conference paper at ICLR 2022
HyperEuler
Buler
Midpoint
RK4
0.5
2
2
2
2 3
2
Figure 8: Displacement variables d and dr of the discretized Timoshenko beam a function of position of the finite elements and time t. The controller optimized with HyperEuler manages to stabilize the beam while the baseline solvers Euler and Midpoint fail, yet requiring less than a third in terms of runtime compared to RK4
Hypersolvers are even more impactful in complex high-dimensional controlled systems.
RELATED WORK
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 & Fatabalian, 2003; Gorissen al., 2010). Differently from these approaches, we investigate a methodology to enable faster solution during 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 explored in (Saveriano et al., 2017) where residual dynamics are learned efficiently along with the control policy while (Fisac 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 residu- als to improve solution accuracy. The synergy between machine learning and optimal control con- tinues 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& & Fatabalian, 2003) pay a computational cost upfront to accelerate downstream sim- ulation. However, ensuring transfer from offline optimization to the online setting is still an open problem. In our approach, we investigate several strategies for 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, approach further lever- ages the latest advances in hardware and machine learning software (Paszke et al., 2019) by solving thousands of ODEs parallel on graphics processing units (GPUs).
CONCLUSION
We presented a novel method for obtaining fast and accurate control policies. Hypersolver mod- els 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 envi- ronments by efficiently solving high-dimensional space-continuous problems.

---

Published as a conference paper at ICLR 2022
CODE OF ETHICS
We acknowledge that all the authors of this work have read and commit to adhering to the ICLR Code of Ethics.
REPRODUCIBILITY STATEMENT
We share the code used in this paper and make it publicly available on Github'. The following ap- pendix 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 Appendix B while in Appendix C we provide details on on the system dynamics, architectures and other experimental details. Additional explanations are also provided as comments in the shared code implementation.
REFERENCES
Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E Rognes, and Garth N Wells. The femics project version 1.5.
Archive of Numerical Software, 3(100), 2015.
Mato Baotić, Francesco Borrelli, Alberto Bemporad, and Manfred Morari. Efficient on-line compu- tation of constrained optimal control. SIAM Journal on Control and Optimization, 47(5):2470- 2489, 2008
Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba. Openai gym, 2016.
John Butcher.
Numerical methods for differential equations and applications. wwt.:atwwww.uuu......eeeeeeerene
22, 12 1997.
Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary dif-
ferential equations, 2019.
S Reynold Chu, Rahmat Shoureshi, and Manoel Tenorio. Neural networks for system identification. IEEE Control systems magazine, 10(3):31-35, 1990.
R. Dormand and P. Prince. A family of embedded runge-kutta formulae. Journal of Computa- tional and Applied Mathematics, 6:19-26, 1980.
Jaime F. Fisac, Anayo K. Akametalu, Melanie N. Zeilinger, Shahab Kaynama, Jeremy Gillula, and Claire Tomlin. general safety framework for learning- based control in uncertain robotic systems, 2018.
Räzvan Florian. Correct equations for the dynamics of the cart-pole system. 08 2005.
C. Garcia, D. M. Prett, and M. Morari. Model predictive control: Theory and practice a survey. Autom,, 25:335-345 1989.
Dirk Gorissen, Ivo Couckuyt, Piet Demeester, Tom Dhaene, and Karel Crombecq. A surrogate modeling and adaptive sampling toolbox for computer based design. The Journal of Machine Learning Research, 11: 11:2051-2055, 2010.
Radek Grzeszczuk, Demetri Terzopoulos, and Geoffrey Hinton. Neuroanimator: Fast neural net- work emulation and control of physics-based models. Proceedings of the 25th annual confer- ence on Computer graphics and interactive techniques, pp. 9-20, 1998.
K.J. Hunt, D. Sbarbaro, R. Žbikowski, and P.J. Gawthrop. Neural networks for control sys- tems-a survey. Automatica, 108-11111 1992. ISSN 0005-1098. doi: https://doi.org 00:001660005-1008309999905532 URL https: httpa arricle/p11/00051082220000
/ www.wciieneei.ee..........m...ee com/ comm/science/
'Supporting reproducibility code is
reaarr//zzee/aaatr//pprsssddeee control
10

---

Published as a conference paper at ICLR 2022
Doug L James and Kayvon Fatahalian. Precomputing interactive dynamic deformable scenes. ACM Transactions on Graphics (TOG), 22(3):879-888 -87888977888888 2003.
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization, 2017.
Frank Lewis, Draguna Vrabie, and Vassilis L Syrmos. Optimal control. John Wiley & Sons, 2012.
Chin-Teng Lin, C. S. George Lee, al. Neural-nettorr-basee fuzzy logic control and decision system. IEEE Transactions on computers, 40(12):1320-1336 1991.
Yinghua Lin and George A Cunningham. A new approach fuzzy-neural system modeling. IEEE Transactions on Fuzzy systems, 3(2):190-198) 1995.
Alessandro Macchelli and Claudio Melchiorri. Modeling and control of the timoshenko beam, the distributed port hamiltonian approach. SIAM J. Control. Optim,, 43:743-767, 2004.
Stefano Massaroli, Michael Poli, Sho Sonoda, Taji Suzuki, Jinkyoo Park, Atsushi Yamashita, and
Hajime Asama. Differentiable multiple shooting layers. CoRR, Abs/2106.0388, 2021. URL fitps::/aaxet .ttriiiiorrg//ooo wtt/abss////ooooooooo.o..o....sr /2106 03885.
David Mayne and Hannah Michalska. Receding horizon control of nonlinear systems. In Pro- ceedings of the 27th IEEE Conference on Decision and Control, pp. 464-465 IEEE, 1988.
Jacopo Panerati, Hehui Zheng, Siqi Zhou, James Xu, Amanda Prorok, and Angela P. Schoellig. Learning to fly a gym environment with pybullet physics for reinforcement learning of multi- agent quadcopter control. CoRR, abs/2103.02142, 2021. URL httts:///aaiiiiii wwtp abs/ 211 2104.0212..
02142.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high- performance deep learning library. arXiv preprint arXiv: 1912.0/72...3 2019.
Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Hyper- solvers: Toward fast continuous-depth models. In H. Larochelle, M. Ranzato, Hadsell, M. F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 21105-21117, Curran Associates, Inc., 2020a. URL https : /ptpe/eeinos....oeo/eeeppo....se neurips. cc/paper/2020/ENa htt.a//eee.166eooae2333333333330 7abbbb Bar000088P-ppr.... pdf.
Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Torchdyn:
A neural differential equations library, 2020b.
Radoslaw Pytlak. Numerical methods for optimal control problems with state constraints. Springer, 2006.
Anil V Rao. A survey of numerical methods for optimal control. Advances in the Astronautical Sciences, 135()))477-528 2009.
Christian Robert and George Casella. Monte Carlo statistical methods. Springer Science & Business Media, 2013.
I Michael Ross and Fariba Fahroo. Issues in the real-time computation of optimal control. Mathe- matical and computer modelling, 43(9-10): 11211 188, 2006.
Matteo Saveriano, Yuchao Yin, Pietro Falco, and Dongheui Lee. Data-efficient control policy search using residual dynamics learning. In 20/7 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4709-47155 2017. doi: 00.1109/ROS.2017.20 2026.33.
Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wet- zstein. Implicit neural representations with periodic activation functions, 2020.
Andrew Taylor, Andrew Singletary, Yisong Yue, and Aaron Ames. Learning for safety-critical control with control barrier functions, 2019.

---

Published as a conference paper at ICLR 2022
Ch. Tsitouras. Runge kutta pairs of order 5(4) satisfying only the first column simplifying assump- Kunte-kuttaaa tion. Computers Mathematics with Applications, 62(2):770-775, 775, 2011. ISSN 0898-12881 1221. doi: https://bbi.org//..1106.... (0006ff.cmmaa22011066600000 URL htts:: / www. sciencedirect. com/
sci / 1004766.
S Vadali, Hanspeter Schaub, and K Alfriend. Initial conditions and fuel-optimal control for forma- tion flying of satellites. In Guidance, Navigation, and Control Conference and Exhibit, pp. 4265, 1999.
Yue J Zhang, Andreas A Malikopoulos, and Christos G Cassandras. Optimal control and coordina- tion of connected and automated vehicles at urban traffic intersections. In 2016 American Control Conference (ACC), pp. 6227-62222222222 6232.. IEEE, 2016.
Liu Ziyin, Tilman Hartwig, and Masahito Ueda. Neural networks fail to learn periodic functions and how to fix it, 2020.
12

---

Published as a conference paper at ICLR 2022
A ADDITIONAL HYPERSOLVER MATERIAL
A.1 EXPLICIT HYPERÉULER FORMULATION
Our analysis in the experiments takes into consideration the hypersolved version of the Euler scheme, namely HyperEuler. Since Euler is first-order method, it requires the least number of function evaluations (NFE) of the vector field f in (1) and yields second order local trun- cation error éx == Rx This error larger than other fixed-step solvers and thus has the most room for potential improvements. The base solver scheme th of (4) can be written as bb (42, Th, us) = f (1a, :R, ux), which is approximating the next state by adding an evaluation of the vector field multiplied by the step size e. We can write the HyperEuler update explicitly as
FR+1 = 2 + ef Zh, un) + e² gw (t&, DRy uk))
(12)
while we write its residual as
1 te, - x(tx) ef 2 uk))
R (x((x), ((t)))) 2
(13)
A.2 HYPERSOLVERS FOR TIME-INVARIANY SYSTEMS
time-invariant system with invarinvariant controller can be described as following
((t) = f (x((t, u(x(t)))
(14)
x00 3op
in which f and u do not explicitly depend on time. The models considered in the experiments satisfy this property.
B CONTROL POLICY DETAILS
B. OPTIMAL CONTROL COST FUNCTION
The general form of the integral cost fanctional can be written as follows
J(x(t), u(t)) = [z T (ty) (aa))=x"0 x²|p[=[t,) - ==1+ /'' (Ix (I==())=)Q))z - -r'] + u (t))t(t) dt
(15) where matrix P is a penalty deviations from the target x² of the last states, penalizes all deviations from the target of intermediate states while R is a regulator for the control inputs. Eval- uation of (15) usually requires numerical solvers such the proposed hypersolvers of this work. Discretizations of the cost functional are also called cost finction in the literature.
B.2 MODEL PREDICTIVE CONTROL FORMULATION
The following problem is solved online and iteratively until the end of the time span
T-1
min
J (Ik, ux))
k=0
subject to
(() =
u(())
(16)
x0) 200
tET
where J is cost function and T is the receding horizon over which predicted future trajectories are optimized.
B.3 NEURAL OPTIMAL CONTROL
We parametrize the control policy of problem (2) as ug free parameters. Specifically, we consider the case of neural optimal control in which controller us is multi-layer perceptron. The optimal control task is to minimize the cost function J described in (15) and we do so by optimizing the parameters via SGD; in particular, we use the Adam (Kingma & Ba, 2017) optimizer for all the experiments.
t, (x(t, where 0 is a finite set of
13

---

Published as a conference paper at ICLR 2022
C EXPERIMENTAL DETAILS
In this section we include additional modeling and experimental details divided into the different dynamical systems.
C.I HYPERSOLVER NETWORK ARCHITECTURE
We design the hypersolver networks Iw forwwrd neural networks. Table summarizes the parameters used for the considered controlled systems, where Activation denotes the activation func- tions, i.c. SoftPlus: z + log(1 + e²), Tanh: z + ="
e--e-" and Snake - + sin (ax)
(Ziyin et al., 2020). We also use the vector field f as an input of the hypersolver, which does not
Table I: Hyper-ppaameetrs for the hypersolver networks in the experiments.
Spring-Man Mass Inverted Pendulum
Cart-Pole
Quadcopter
Timoshenko
Beam
Input Layer
5
5
28
322
Hidden Layer
32
32
64
256
Activation
Softplus
Softplus
Snake
Softplus
Snake
Hidden Layer
32
32
32
64
256
Activation 2
Tanh
Tanh 2
Snake
Softplus 12
Snake
Output Layer
2
160
require further evaluation since it is pre-evaluated runtime by the base solver th. We empha- size 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 hy- persolvers also depends on their complexity. Keeping their neural network small enough guarantees that evaluating the hypersolvers is cheaper than resorting more complex numerical routines.
C.2 MPRING-MASS SYSTEM
System Dynamies The spring-mass system considered described in the Hamiltonian formula- tion by
1/m 0
+
(17)
P
where m = 1 [Kg] and & = 0.5 [N/m].
Pre-training methods comparison We select {(x, u) as a uniform distribution with support in X' x X where X' [-20 20] [-20, 20) while 21 = [-100, 100). Nominal solutions 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 different training methods local residual for step size c 0.03 s: stochastic exploration and active error minimization. The optimizer used is Adam with learning rate of 10-3 for 10a epochs.
Hypersolvers with different step sizes We select E(x,, u) as uniform distribution with support in x l where X' = on the accurate system using dopri5 with relative and absolute tolerances set to 10-7 and 10 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 10a epochs.
5-5,5) [-5,
[-5,5] while 4 = 1-20, 20). Nominal solutions are calculated
Thiis architecture refers the optimal control experiment. Details hypersolver models for the general- ization experiment on the inverted pendulum are available Appendix C.3.
In the multi-stage hypersolver experiment consider both the inner stage h and the outer stage Iw with the architecture and jointly trained (more information and ablation study in Appendix C.4).
14

---

Published as a conference paper at ICLR 2022
10
101
Mran Ri-icdual &
Moan Readual
10
10
10
10
10-2
10-3
10
[N]
Euler
Tanh [F3E]
Snake [HC]]
Fuler
[HR] ReLU [HE]
Midpoins
(HR) Training Region
Midpoint
Trainnn Step
SIREN (HF)
RK4
SIREN [HE]
Figure 9: Generalization with different hypersolver activation functions (HyperEuler models marked with "HE") on the inverted pendulum. [Left] Generalization for the controller space outside of the training region (red area). The architecture with Snake can to generalize better compared to other hypersolvers. [Right]
Generalization for different time steps outside of the training step 0.1 a (red line). HyperEuler is able to improve the baseline Euler solver performance even for unseen
C.3 INVERTED PENDULUM
System Dynamics via the following:
We model the inverted pendulum with elastic joint with Hamiltonian dynamics
0 1/m [ -8/m
p
mgl sin
+
(18)
where m 1 [Kg], 0.5 (N/ rad], r [m). 3 0.01 (Ns/ rad]. g 9.81 [m/s²].
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' = [-2T, 2 x 122 2m] and [-10, 10] for the different architectures. We choose as a common time step networks are trained for 100000 epochs with the Adam optimizer and learning rate of 10 The network architectures share the same parameters as the inverted pendulum ones in while the activation functions are substituted by the ones in Figure 4. The SIREN architecture is chosen with hidden layers of size 64. Figure provides an additional empirical results on generalization properties across controller values and step sizes: notice how Snake can generalize to unseen control values better compared to other hypersolvers.
0.1 s; the
Additional visualization Figure provides an additional visualization of the inverted pendulum controlled trajectories from Figure with positions g and momenta over time.
HyperEuler
Euler
Midpoint
RK4
400
2
5
Figure 10: Controlled trajectories of the inverted pendulum with controllers optimized via different solvers.
C.4 CART-POLE
System Dynamics We consider a continuous version a cart- pole system additionally taking into account the full dynamic model in Florian (2005). This formulation considers the friction coefficient
15

---

Published as a conference paper at ICLR 2022
between the track and the cart He inducing force opposing the linear motion well as the friction generated between the cart and the pole Ap, whose generated torque opposes the angular motion. The full cart-pole model is described by the four variables 2, á 0, 0 and the accelerations update is as following
No (mc + mp) 9 mp! (ösin O 0² cos
- 1-- -u 102 (sin t+ sgn( (N.)) cos 0) + Mcg sgg
g sin 8 + cos 0
mc+mp
my.
4
mp oos
my (cos Me sgn(N.
(19)
u + mpl ( sin 0 0 cas 0 Mo 3
m + mp
where m. g]. mp 0.1 [Kg]. 0.5 [m] and g normal force acting on the cart. For simulation purposes, we consider its sign to be always positive when evaluating the sign (sgn) function as the cart should normally not jump off the track. Setting Me. Up to 0 results in the same dynamic model used in the OpenAl Gym (Brockman et al., 2016)
9.81 [m/s²]. N₂ represents the
implementation.
20-1
SheSte MMB
kointtoiiinng MS MS trinnng MS nooiieer
MB Residul dynamics
Edler (nominal)
Milpoon (ronioa)) Midpoin (raccurat)
Figure 11: One step Mean Absolute Frror (MAF) for multi-stage hypersolvers and different solvers as well as correction schemes in the ablation study. Multi-stage Hypersolver (MS) with joint training and Midpoint base solver iterating an inaccurate vector field agnostic of friction forces outperforms the Midpoint solver with full knowledge of the vector field.
Multistage training strategies We study two different training strategies for the inner and outer networks An and f in (11). We first consider a joint training strategy in which both stages are trained the same time via stochastic exploration. Secondly, we do a separated training in which only the inner stage network hn is trained first and then the outer stage network g. is added and only its parameters are trained a finetuning process. We find that, as shown in Figure 11, joint training yields slightly better results. A further advantage of jointly training both stages is that only a single training procedure is required.
Ablation study We consider the same model as Multi Hypersolver with base Midpoint solver but no first-stage fg. which corresponds to learning the residual dynamics only, and we train this model with stochastic exploration. We show in Figure 11 that while the residual dynamics model can improve the one-step error compared to the base solver on the inaccurate dynamics, it performs worse the Multi-stage Hypersolver scheme. We additionally study the contribution of each stage in the prediction error improvements by separately zeroing out the contributions of the inner and outer stage. While iterating over the inner stage only improves on the base-solver error, including the outer stage further contributes in improving the error. We notice how the excluding the inner- stage yields higher errors: this may be due to the fact that the anner-staee specializes in correcting the first-order vector field inaccuracies while the outer-stage corrects the one step base solver residual.
Additional experimental details For the Multi-stage Hypersolver control experiment, we pre- train both inner and outer stage networks ht and gw in (11) at the same time using stochastic exploration. The base solver is chosen as the second-orde namies (19) with Mes Ho set to O. The nominal dynamics considers non-null friction forees: we set the cart friction coefficient to H. = 0.1 and the one of the pole to 14, = 0.03. We note how the friction coefficients make the vector field (19) non-smooth: simulation through adaptive-step size
Midpoint iterating on the partial dy-
16

---

Published as a conference paper at ICLR 2022
solvers tsit5 results experimentally time-consuming, hence resort to RK4 for training the hypersolver networks. Nonetheless, as shown in the error propagation of Figure 12, this does not degrade the performance of the trained multi- stage hypersolver scheme. All neural networks in the experiements, including the ablation study, are trained with the Adam optimizer for x 10 epochs, where we the learning rate 10-2 for the first 10a epochs, then decrease it 10-3 for 10 epochs and to 10 for the last 10a.
C.5 QUADCOPTER
System Dynamics update in (Panerati al., 2021) for batched training in PyTorch. The following accelerations update describes the dynamic model
The quadcopter model is a suitably modified version the explicit dynamics
R [0, 00 Σ -u 10,0, m.))) m
(20)
& = J-1 (T(1, kFF lo, (wo: wp w3 w ab (J
mwhere = [xr y, z] corresponds the drone positions and vb [o, 0, its angular positions; R and J are its rotation and inertial matrices respectively, ---) is a function calculating the torques induced by the motor speeds wi. while arm length I, mass m, gravity acceleration constant along with k and ky are scalar variables describing the quadcopter's physical properties.
Position
Velucity
Angular Position
Angular Velocity
10
10
10-2
SMAPE
SMAPP
SMATE
10-3
10-4
Time 6]
Time
Time
Time
Multinage Hyparachor finaccurato modal)
Midguiet finaccumno moddl)
Edler (suminal mode))
Midpoint (nomma) modal)
Figure 12: Symmetric Mean Absolute Percentage Error (SMAPE) propagation along controlled trajectories the cart-pole system. The Multi-stage Hypersolver knowledge limited to inaccurate model manages to outperform the Euler solver iterating on the accurate dynamics in terms of positions and angular positions.
C.6 TIMOSHENKO BEAM
System Dynamics We consider as a system from the theory of continuum dynamics the Tim- oshenko beam with no dissipation described in (Macchelli & Melchiorri, 2004; Massaroli et al., 2021). The system can be described in the coenergy formulation by the following partial differen- tial equation (PDE)
PA
0
0
0
Vt
0
0
0
0
Up
0
1
0
0
0
Ur
0
0
0
1
Ur
(21)
0
0
C
0
Ot
Tr
0
0
0
0
Or
0
0
C.
a
where p is the mass density. A is the cross section area, 4, is the rotational inertia, C, and C2 are the shear and bending compliance; the discretized state variables Ut, V,- represent the translational and rotational velocities respectively while dz, a, denote the translational and rotational displacements. cantilever beam. order to discretize the problem, we implement software routine based on the fenics (Alnæs et al., 2015) open-source software suite to obtain the finite-llements discretization of the Timoshenko PDE of (21) given the number of elements, physical parameters of the model and initial conditions of the beam. We choose 40 elements discretization of the PDE for total of 160 dimensions the discretized state = [0r, Uy. ot, o,] and we initialize the beam at time t = =(x,00
sin( (3Tx), 0, 0). The system thus be reduced to the following controlled
17

---

Published as a conference paper at ICLR 2022
linear system
-MAD PA -MI,
Vt
X
M BF
Up
X
x
MM
Do
M
Br
+
dt
Mc. D2
x
6r_ LM Di
Mc. Do
Or
x
(22)
where the mass matrices MpA: MI., Mc.: Mc,. matrices Do: D1, D2. vectors Bp, BT are computed through the fenics routine and boundary controllers ug and uy are the control torque and the control force applied 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' x U. We restrict the boundary control input values = [ug. 2 in [-1, x [-1,1]. As for the state space, naively generating probability distribution with box boundaries on each of the 160 dimensions of X' would require an inefficient every combination is physically feasible due to the Timoshenko beam's s structure. We solve this problem by propagating batched trajectories with RK4 from the initial boundary condition 2(x, 0)
search over this high-dimensional space: in fact,
sin ,, 0, with random control actions sampled from uniform distribution with support in [-1, 1] [-1, applied for time I ((0.002,1) 1] 8. We save the states (=(x, (1) and forward propagate from these states again by sampling from the controller and time distributions. We repeat the process K times and obtain sequence [x t1 ) {a(x, {x]']] ]' of batched initial conditions characterized by physical fea- sibility. Finally, train the hypersolver with stochastic exploration by sampling from the generated distribution {xx, u) on local one-step residuals as deseribed in Section 5.4. This initial state gener- ation strategy is repeated every 100 epochs for guaranteeing 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.
10
10a
10a
10-3
[(6]
13
[16]
Euler
HyperEaler
Midpoint
R344
Figure 13: Mean Absolute Error (MAE) propagation on velocities Vt, e, and displacements o,, for the finite elements of the discretized Timoshenko beam along controlled trajectories. While solutions from Euler and Midpoint quickly diverge due the system's 's stiffness, HyperEuler manages only contain errors but even outperform the fourth-order RK4 whilst requiring fraction of number of vector field evaluations.
Additional details the results We report additional details regarding the runtime of the exper- iments on the Timoshenko Beam. As for the training time, training the hypersolver for 10 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) training iteration, Euler (2.01 + 0.04) Midpoint (4.02 + 0.08) and RK4 (8.24 = 0.14) s. Experiments were run on the CPU of the machine described in Section C.7.
C.7 HARDWARE AND SOFTWARE
Experiments carried out machine equipped with AMD RYZEN THREADRIPPER 3960X CPU with 48 threads and two NVIDIA RTX 3090 graphic cards. Software-win we used PyTorch (Paszke al., 2019) for deep leamiing and the torchdyn (Poli et al., 2020b) and
18

---

Published as a conference paper at ICLR 2022
torchdiffeq (Chen et al., 2019) libraries for ODE solvers. We additionally share the code used in this paper and make it publicly available on Github".
'Supporting reproducibilltyy code is at
httpa //github.co....to//e/ www.DDiffEgM../afffee.
reearca/zzeenntteer/ppppaaaeeeee control
19


Reviews:
[{"review_id": "RY2aR8fPCal", "reviewer": "Reviewer_RBHY", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths: I thought the paper was a crisp and clear read, and that all technical components (both novel and tutorial) and all experiments were both sound and very well-explained.\n\nWeaknesses: The systems that are tested on are relatively small scale/low-dimensional. Given that hypersolver frameworks are presumably most necessary in larger scale systems that are indeed expensive to solve, this is in my view a large limitation in terms of understanding the value proposition of the present work. It would greatly strengthen the present work if demonstrations on such large-scale/high-dimensional systems could be shown.\n\nMinor points (not affecting my score):\n* $\\psi_\\epsilon$ in Equation (3) could be more explicitly defined, as it is of course a critical variable for the rest of the paper.\n* The differentiation of the present work from the prior hypersolvers work could be made more clear, and in general the writing could better signpost what aspects of the present work are \"new\" contributions. If I understand correctly, the major distinction from the prior work is that the prior work considered uncontrollable systems, whereas the present work considers controllable systems, which yields some differences both in terms of how the hypersolver is trained and in how it is used. However, this should be spelled out more explicitly.", "summary_of_the_paper": "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.", "summary_of_the_review": "I thought this paper was a sound, enjoyable, and clear extension of the previous work on hypersolvers to the setting of controlled dynamical systems. However, providing experiments on larger-dimensional systems would in my view significantly strengthen the paper by more fully demonstrating the potential benefits of the proposed method. That said, I am marginally inclined to accept the paper nonetheless rather than rejecting it solely on that basis.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "Jl7-q5678ku", "reviewer": "Reviewer_J9x8", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I thought this paper on the whole was nicely written and easy to understand. \n\nMy main issue with the paper is the limited application and its relevance to the ICLR community. \nThe idea of using NNs to warm start ODE solvers is relevant, but this idea is presented in the paper from Poli et al from Neurips 2020 (althought the paper refence of this work is missing the venue).\nHowever the specific application to optimal control ODEs is somewhat limited and would probably be better placed at a control conference like CDC, ACC or L4DC. \nThis is reinforced by the experimental settings being quite simple and lwo dimensional. \nWhile ODE-based optimal control is an important setting, the ML community typically use other methods for continuous control. \n\nI also had some minor issues with mathematical clarity. For example the losses in Eq 6 and 6 dont have inputs. In the (unnumbered) equation in 4.2 and Eq. 8 , it is not clear where w enters the objective.  ", "summary_of_the_paper": "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. \nThis warm-starting improves downstream accuracy when the computational budget is fixed. \nThis method is applied for trajectory optimziation for pendulum swing-up and MPC for cartpole and quadcopter systems.", "summary_of_the_review": "The application of hypersolves to optimal control is somewhat incremental and not of great relevance to the ICLR community. I woudl be willing to change my score if the reviewers and AC disagree. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "DQ118GRM-7", "reviewer": "Reviewer_jwFx", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper proposes to apply the idea of hypersolvers to numerical optimal control. Through learning the truncation residual dynamics, the model is able to approximate the actual dynamics with a combination of low-order ODE solver with learned truncation residual dynamics. Further, if the system model is not perfectly accurate, a multi-stage hypersolver can be used to separately learn the system residual dynamics and truncation residual dynamics. In experiments, the method is shown to perform on par with the higher-order ODE solvers. Also, experiments show that certain activation functions (SIREN and Snake) perform better than the standard ones (Tanh and ReLU).\n\nHere are some of my detailed comments:\n1. In (3), \\epsilon appeared before it is introduced in (4), hence making it undefined and unclear.\n2. In section 3.2, maybe more explanation can be given on how (2) is solved?\n3. In the pendulum experiment, the hypersolver approach achieves similar performance in terms of final positions when compared with Midpoint and RK4, but using similar number or 50% less of FLOPS. For the other two experiments, I don’t seem to find how does the method compares in terms of FLOPS or runtime. As in the pendulum case, the proposed hypersolver approach does not outperform a lot in terms of runtime or FLOPS.\n4. In the pre-training of the hypersolvers, it seems to require a very large number of epochs, maybe more information can be provided on how many trajectories are used for training and how long the training time is.\n5.  In section 5.2, how does the tolerance of the ODE solver affect the experiment results, maybe some ablation studies and explanations can be included to help understand.\n6. For the multi-stage hypersolvers, there are some literature on learning the system residual dynamics for control, such as [1][2][3]. Also some ablation studies can be done by comparing the multi-stage hypersolver approach with just learning system residual dyanamics.\n\n[1] Saveriano, Matteo, Yuchao Yin, Pietro Falco, and Dongheui Lee. \"Data-efficient control policy search using residual dynamics learning.\" In 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 4709-4715. IEEE, 2017\n\n[2] J. F. Fisac, A. K. Akametalu, M. N. Zeilinger, S. Kaynama, J. Gillula, and C. J. Tomlin. A general safety framework forlearning-based control in uncertain robotic systems.IEEE Transactions on Automatic Control, 64(7):2737–2752,2018.\n\n[3] A. Taylor, A. Singletary, Y. Yue, and A. Ames. Learning for safety-critical control with control barrier functions. InLearning for Dynamics and Control, pages 708–717. PMLR, 2020.\n", "summary_of_the_paper": "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.", "summary_of_the_review": "Overall, the paper proposes to use hypersolvers to speed up ODE solving for numerical control. The idea is promising but it needs more experiment results comparisons, especially in terms of runtime, to justify the efficacy of using hypersolvers for optimal control. Also, the paper could consider expanding the explanation on how the optimal control problem (2) is solved in their experiments. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors propose a novel hypersolver framework for solving numerical optimal control problems, learning a low order ODE and a neural network based residual dynamics. They compare their framework with traditional optimal control solvers on a number of control tasks and demonstrate superior performance.

The reviewers are in consensus that the paper makes significant contributions that are validated by the experimental results. The only concern was that the experiments are largely on low dimensional systems, but the reviewers agreed that the results are still worthy of acceptance.
Decision: Accept (Poster)