Learning to Control PDEs with Differentiable Physics
Authors: Philipp Holl
Nils Thuerey
Vladlen Koltun
Venue:  (2020)
Keywords: Differentiable physics
Optimal control
Deep learning

Abstract: Predicting outcomes and planning interactions with the physical world are long-standing goals for machine learning. A variety of such tasks involves continuous physical systems, which can be described by partial differential equations (PDEs) with many degrees of freedom. Existing methods that aim to control the dynamics of such systems are typically limited to relatively short time frames or a small number of interaction parameters. We present a novel hierarchical predictor-corrector scheme which enables neural networks to learn to understand and control complex nonlinear physical systems over long time frames. We propose to split the problem into two distinct tasks: planning and control. To this end, we introduce a predictor network that plans optimal trajectories and a control network that infers the corresponding control parameters. Both stages are trained end-to-end using a differentiable PDE solver. We demonstrate that our method successfully develops an understanding of complex physical systems and learns to control them for tasks involving PDEs such as the incompressible Navier-Stokes equations.

Published as a conference paper at ICLR 2020
LEARNING TO CONTROL PDEs WITH DIFFERENTIABLE PHYSICS
Philipp Holl
Vladlen Koltun Intel Labs
Nils Thuerey
Technical University of Munich
Technical University of Munich
ABSTRACT
Predicting outcomes and planning interactions with the physical world are long- standing goals for machine learning. A variety of such tasks involves contin- uous physical systems, which can be described by partial differential equations (PDEs) with many degrees of freedom. Existing methods that aim to control the dynamics of such systems are typically limited to relatively short time frames or a small number of interaction parameters. We present novel hierarchical predictor- corrector scheme which enables neural networks learn to understand and control complex nonlinear physical systems over long time frames. We propose to split the problem into two distinct tasks: planning and control. To this end, we intro- duce predictor network that plans optimal trajectories and control network that infers the corresponding control parameters. Both stages are trained end-to-end using differentiable PDE solver. We demonstrate that method successfully develops an understanding of complex physical systems and learns to control them for tasks involving PDEs such as the incompressible Navier-Stokes equations.
INTRODUCTION
Intelligent systems that operate in the physical world must be able to perceive, predict, and interact with physical phenomena (Battaglia al., 2013). In this work, we consider physical systems that can be characterized by partial differential equations (PDEs). PDEs constitute the most fundamen- tal description of evolving systems and are used to describe every physical theory, from quantum mechanics and general relativity to turbulent flows (Courant & Hilbert, 1962; Smith, 1985). We aim to endow artificial intelligent agents with the ability to direct the evolution of such systems continuous controls.
Such optimal control problems have typically been addressed via iterative optimization. Differ- entiable solvers and the adjoint method enable efficient optimization of high-dimensional sys- tems (Toussaint et al., 2018; de Avila Belbute-Peres et al., 2018; Schenck & Fox, 2018). However, direct optimization through gradient descent (single shooting) test time is essource-intensiv and may be difficult to deploy in interactive settings. More advanced methods exist, such as multiple shooting and collocation, but they commonly rely modeling assumptions that limit their applica- bility, and still require computationally intensive iterative optimization at test time.
Iterative optimization methods expensive because they have to start optimizing from scratch and typically require a large number of iterations to reach an optimum. In many real-world control prob- lems, however, agents have to repeatedly make decisions in specialized environments, and reaction times are limited to a fraction of second. This motivates the use of data-driven models such as deep neural networks, which combine short inference times with the capacity build an internal representation of the environment.
We present a novel deep learning approach that can learn to represent solution manifolds for a given physical environment, and is orders of magnitude faster than iterative optimization techniques. The core of our method is hierarchical predictor-correate scheme that temporally divides the problem into easier subproblems. This enables us to combine models specialized to different time scales in order to control long sequences of complex nigh-dimensional systems. We train our models using a differentiable PDE solver that can provide the agent with feedback of how interactions at any point in time affect the outcome. Our models learn to represent manifolds containing a large number of solutions, and can thereby avoid local minima that can trap classic optimization techniques.

---

Published as a conference paper at ICLR 2020
We evaluate method on variety of control tasks in systems governed by advection- -iffiioion PDEs such as the Navier-Stokes equations. We quantitatively evaluate the resulting sequences on how well they approximate the target state and how much force was exerted on the physical system. Our method yields stable control for significantly longer time spans than alternative approaches.
2 BACKGROUND
Physical problems commonly involve nonlinear PDEs, often with many degrees of freedom. this context, several works have proposed methods for improving the solution of PDE problems (Long et al., 2018; Bar-Sinai et al., 2019; Hsich et al., 2019) or used PDE formulations for unsupervised optimization (Raissi et al., 2018). Lagrangian fluid simulation has been tackled with regression forests (Ladicky et al., 2015), graph neural networks (Mrowca et al., 2018; Li et al., 2019). and con- tinuous convolutions (Ummenhofer al., 2020). Data-driven turbulence models were trained with MLPs (Ling et al., 2016). Fully-convolutional networks were trained for pressure inference (Tomp- son al., 2017) and advection components were used in adversarial settings (Xie al., 2018). Temporal updates in reduced spaces were learned via the Koopman operator (Morton et al., 2018). In a related area, deep networks have been used to predict chemical properties and the outcome of chemical reactions (Gilmer et al., 2017; Bradshaw et al., 2019).
Differentiable solvers have been shown to be useful in variety of settings. Degrave et al. (2019) and de Avila Belbute-Pered et al. (2018) developed differentiable simulators for rigid body me- chanics. (See Popovic et al. (2000) for earlier work in computer graphics.) Toussaint et al. (2018) applied related techniques to manipulation planning. Specialized solvers were developed to infer protein structures (Ingraham et al., 2019), interact with liquids (Schenck & Fox, 2018), control soft robots (Hu et al., 2019), and solve inverse problems that involve cloth (Liang et al., 2019). Like ours, these works typically leverage the automatic differentiation of deep learning pipelines (Griewank & Walther, 2008; Maclaurin al., 2015; Amos & Kolter, 2017; Mensch & Blondel, 2018; van Merriënboer et al., 20:8; Chen et al., 2018; Bradbury et al., 2018; Paszke et al., 2019; Tokui et al., 2019). However, while the works above target Lagrangian solvers, i.e. reference frames moving with the simulated material, we address grid-based solvers, which are particularly appropriate for dense, volumetric phenomena.
The adjoint method (Lions, 1971; Pironneau, 1974; Jameson, 1988; Giles & Pierce. 2000; Bewley, 2001; McNamara et al., 2004) is used by most machine learning frameworks, where it is commonly known as reverse mode differentiation (Werbos, 2006; Chen et al., 2018). While a variety of spe- cialized adjoint solvers exist (Griewank et al., 1996; Fournier et al., 2012; Farrell et al., 2013), these packages do not interface with production machine learning frameworks. A supporting contribu- tion of our work a differentiable PDE solver called Frow that integrates with TensorFlow (Abadi et al., 2016) and PyTorcc (Paszke et al., 2019). It is publicly available at https://gituub..om/um pbs/PhiFlow.
3 PROBLEM
Consider physical system u(x 1) whose natural evolution is described by the PDE
du
P ( Ou Gd
dx Ox² 2....., v(t))
Ot
where P models the physical behavior of the system and (t) denotes external factors that can influence the system. We now introduce an agent that can interact with the system by controlling certain parameters of the dynamics. This could be the rotation of motor or fine-grained control over field. We factor out this influence into force term F, yielding
au = P u. Øt
Tu D2u
+ F(!).
(2)
Oz dx
The agent can now be modelled as a function that computes F(t). As solutions of nonlinear PDEs were shown to yield low- dineennniional manifolds (Foias et al., 1988; Titi, 1990), we target solution manifolds of Ft) for a given choice of P with suitable boundary conditions. This motivates our choice to employ deep networks for agents.
2

---

Published as a conference paper at ICLR 2020
In most real- world scenarios, it is not possible observe the full state of physical system. When considering a cloud of smoke, for example, the smoke density may observable while the velocity field may not be seen directly. We model the imperfect information by defining the observable state of u as o(u). The observable state problem dependent, and our agent is conditioned only on these observations, i.e. it does not have access the full state
Using the above notation, we define the control task as follows. An initial observable state og of the PDE well a target o" are given (Figure la). We are interested in reconstructed trajectory (()) that matches these states at to and l.,, i.e. On ((u()))) o* o(u(t.)). and minimizes the amount of force applied within the simulation domain D (Figure lb):
fm S. [F)( dx dt.
1F(a((]=
Taking discrete time steps At, the reconstructed trajectory u a sequence of n = (f. - 4/) states.
(a) Task
(b) Trajectories
Observation
Observation
When an observable dimension cannot be con- trolled directly, there may not exist any trajec- tory utt) that matches both Op and o". This stem from either physical constraints merical limitations. In these cases, we settle for an approximation of o*. To measure the quality of the approximation of the target, we define an observation loss Lg The form of this loss can be chosen to fit the problem. We combine Eq. 3 and the observation loss into the objective func- tion
Op
Hidden state h
Figure I: Illustration of possible trajectories. The grey lines represent the unperturbed evolution of the physical system. The amount of applied force corresponds to how far the trajectory deviates from the natural evolution.
L[u((t) = a. LF(u(()) + (t.()t))) (4)
with > O. We square brackets to denote functionals, i.e. functions depending fields or series rather than single values.
PRELIMINARIES
Differentiable solvers. Let (ax 1) be described by PDE in Eq. A regular solver can move the system forward in time via Euler steps:
w(t++1)= Solver(u(t)), y(t+) = u(;) + At ((((t))..,.t()))).
(5)
Each step moves the system forward by time increment At. Repeated execution produces tra- jectory u(t) that approximates a solution to the PDE. This functionality for time advancement by itself is not well-suited to solve optimization problems, since gradients by finite differencing. For high-dimensional or continuous systems, this method becomes computa- tionally expensive because full trajectory needs be computed for each optimizable parameter.
only be approximated
Differentiable solvers resolve this issue by solving the adjoint problem (Pontryagin, 1962) via ana- lytic derivatives. The adjoint problem computes the same mathematical expressions while working vectors. A differentiable solver can efficiently compute the derivatives with
with nower-dimensionss respect to any of its inputs, i.e. @u(t++1)/ 1/0ut(t))t)))) (u.(tt) and du(t++1) (u(tt+1)/dy(t))) This allows for gradient- based optimization of inputs or control parameters over an arbitrary number of time steps.
Iterative trajectory optimization. Many techniques exist that try find optimal trajectories by starting with an initial guess for Ftt) and slightly changing it until reaching an optimum. The sim- plest of these known as single shooting. then backpropagates the loss through the whole sequence optimize the controls (Kraft, 1985; Leineweber et al., 2003). Replacing F(t) with an agent F(tto,,, o"), which can be parameterized by deep network, yields simple training method. For sequence of frames, this setup contains linked copics of the agent and depicted in Figure 2. We refer to such an agent as a control force estimator (CFE).
one optimization step, it simulates the full dynamics,
3

---

Published as a conference paper at ICLR 2020
Optimizing such chain of CFEs is both computationnlly expensive and causes gradients pass through a potentially long sequence of highly nonlinear simulation steps. When the reconstruction is close optimal trajectory, this is not a problem because the gradients Au are small and the operations executed by the solver are differentiable by construction. The solver can therefore be locally approximated by a first-order polynomial and the gradients can be safely backpropagated. For large Au, c.g. at the beginning of an optimization, this approximation breaks down, causing the gradients become unstable while passing through the chain. This instability in the training process can prevent single-shooting approaches from converging and deep networks from learning unless they are initialized near an optimum.
Alternatives to single shooting exist, promising better and more efficient convergence. Multi- ple shooting (Bock & Plitt, 1984) splits the trajectory into segments with additional defect con- straints. Depending on the physical system, this method may have to be adjusted for specific prob- lems (Treuille al., 2003). Collocation schemes (Hargraves & Paris, 1987) model trajectories with splines. While this works well for particle trajectories, it is poorly suited for Eulerian solvers where the evolution of individual points does not reflect the overall motion. Model reduction can be used to reduce the dimensionality nonlinearity of the problem, but generally requires knowledge. When applicable, these methods can converge faster or in a more stable manner than single shooting. However, as we are focusing on general optimization scheme in this work, will use single shooting and its variants as baseline comparisons.
Supervised and differentiable physies losses. One of the key ingredients in training machine learning model is the choice of loss function. For many tasks, supervised losses are used, i.e. losses that directly compare the output of the model for a specific input with the desired ground truth. While supervised losses can be employed for trajectory optimization, far better loss functions are possible when a differentiable solver is available. We will refer to these as differentiable physics loss functions. In this work, we employ combination of supervised and differentiable physics losses, as both come with advantages and disadvantages.
One key limitation of supervised losses is that they can only measure the error of a single time step- Therefore, an agent cannot get any measure of how its output would influence future time steps. Another problem arises from the form of supervised training data which comprises input-output pairs, which may be obtained directly from data generation or through iterative optimization. Since optimal control problems one input. This ambiguity in the supervised training process will lead to suboptimal predictions as the network will try to find compromise between all possible outputs instead of picking one of them.
generally not unimodal, there can exist multiple possible outputs for
Differentiable physics losses solve these problems by allowing the agent to be directly optimized for the desired objective (Eq. 4). Unlike supervised losses, differentiable physics losses require a differentiable solver backpropagate the gradients through the simulation. Multiple time steps can be chained together, which is key requirement since the objective (Eq. 4) explicitly depends on all time steps through L=[u(t)] (Eq. 3). As with iterative solvers, one optimization step for a sequence of n frames then invokes the agent n times before computing the loss, each invocation followed by a solver step. The employed differentiable solver backpropagates the gradients through the whole sequence, which gives the model feedback on (i) how its decisions change the future trajectory and
(a) Forward pass
CFE
Solver
CFE
Solver
CFE
Solver
o(u')
(b) Backward pass
L
O.
CFE
Adjoint
CFF
Adjoint
CFE"
Adjoint
Aut
(c) Weight update
Aw
Aw
Aa
t1
Ex-1t
Cn
Figure 2: Single-shooting optimization with control force estimator (CFE). (a) The forward pass simulates the full sequence. (b) Backpropagation computes the adjoint problem. (c) Weight updates are accumulated and applied to the CFE.
4

---

Published as a conference paper at ICLR 2020
(ii) how handle states imput that were reached because of its previous decisions. Since ground truth needs to be provided, multi-modal problems naturally converge towards one solution.
METHOD
In order to optimally interact with physical system, an agent has to (i) build an internal repre- sentation of an optimal observable trajectory o(u(t)) and (ii) learn what actions to take to move the system along the desired trajectory. These two steps strongly resemble the predictor-correett method (Press et al., 2007). Given o(l), a predictor-corrector method computes o(I + At) in two steps. First, a prediction step approximates the next state, yielding o (1++ At). Then, the correction uses (t + A!) to refine the initial approximation and obtain ol + At). Each step can, to some degree, be learned independently.
This motivates splitting the agent into two neural networks: an observation predictor (OP) network that infers intermediate states (l,),, 21, 1}. planning out trajectory, and corrector network (CFE) that estimates the control force F(t:lo(us o++1) to follow that trajectory as closely as possible. This splitting has the added benefit of exposing the planned trajectory, which would otherwise be inaccessible. As we will demonstrate, it is crucial for the prediction stage to incorporate knowledge about longer time spans. We address this by modelling the prediction as a temporally hierarchical process, recursively dividing the problem into smaller subproblems.
To achieve this, we let the OP not directly infer o (11+1 (((:), , o²) but instead model it to predict the optimal center point between two states at times L; lj with i, E 21,2,...z 1}.j > i, i.e.
(tt + tty //2 o 0,). This function is much more general than predicting the state the next time step since two arbitrary states be passed as arguments. Recursive OP evaluations can then partition the sequence until a prediction of for every time step t; has been made.
This scheme naturally enables scaling to arbitrary time frames or arbitrary temporal resolutions, assuming that the OP can correctly anticipate the physical behavior. Since physical systems often exhibit different behaviors arbitrary time spans, we condition the OP on the time scale it is evaluated on by instantiating and training unique version of the model for every time scale. This simplifies training and does not significantly increase the model complexity as we use factors of two for the time scales, and hence the number of required models scales with Ologg n).. We will refer one instance of an OP, by the time span between its input states, measured in the number of frames = (4g 4)/
different time scales and the OP can be called with states separated by
Execution order. With the CFE and OP, as building blocks, many algorithms for solving the control problem, i.e. for computing F(t), can be assembled and trained. We compared a variety of algorithms and found that scheme we will refer to prediction refinement produces the best re- sults. It is based the following prin- ciples: (i) always use the finest scale OP possible to make a prediction, (ii) execute the CFE followed by a solver step as soon as possible, (iii) refine predictions after the solver has computed the next state. The al- gorithm that realizes these goals shown in Appendix B with an example for n = 8. To understand the algorithm and resulting execution orders, it is helpful to consider simpler algorithms first.
(a) Prediction first
(b) Staggered
(c) Refined
Timet
Time
Time
OP
OP
OP
OP
OP
OP
OP
OP
CFF
CFF
CFFF
CFE
OP
OP
Solver
CFF
Solvet
CFF
Solver
Figure 3: Overview of the different execution schemes.
The simplest combination of CFE and OP,, invocations that solves the full trajectory, shown in Fig- ure 3a, can be described as follows. Initially, all intermediate states are predicted hierarchically. The first prediction is the half-way point o² (l On, o²), generated by the OP,. Using that as input to an OP,// results in new predictions at fn/4; tan/4- Continuing with this scheme, a prediction can be made for each t; i I,..., 1. Next. the actual trajectory is evaluated step by step. For each step tú, the CFE computes the control force Ftt,) conditioned on the state at t; and the prediction o? (1++1)). Once F(..) is known, the solver step the simulation to the next state at 44+1- This al-
5

---

Published as a conference paper at ICLR 2020
gorithm finds trajectory in time (nn) since n CFE calls and n- OP calls required in total (see Appendix B). However, there are inherent problems with this algorithm. The physical constraints of the PDE and potential approximation errors of the CFE can result in observations that only matched partially. This can result in the reconstructed trajectory exhibiting undesirable oscillations, often visible jittering. When subsequent predictions do not line up perfectly, large forces may be applied by the CFE or the reconstructed trajectory might stop following the predictions altogether.
This problem can be alleviated by changing the execution order of the two-stage algorithm described above. The resulting algorithm is shown in Figure 3b and will be referred to as staggered execution. In this setup, the simulation advanced as soon as a prediction for the next observable state exists and OPs are only executed when their state at time t; is available. This staggered execution scheme allows future predictions to take deviations from the predicted trajectory into account, preventing a divergence of the actual evolution o(u(()) from the prediction o² (t).
While the staggered execution allows most predictions correct for deviations from the predicted trajectory o", this scheme leaves several predictions unmodified. Most notably, the prediction (4n72), which is inferred from just the initial state and the desired target, remains unchanged. This prediction must therefore be able to guide the reconstruction in the right direction without knowing about deviations in the system that occurred up to /2/2-1. As a practical consequence, a network trained with this scheme typically learns to average over the deviations, resulting in blurred predictions (see Appendix D.2).
Algorithm 1: Recursive algorithm computing the prediction refinement. The algorithm is called via Reconstruct Ow, absent] to reconstruct a full trajectory from On to o..
function Reconstruct (o( un); on O2n]:
Input Initial observation of Un), observation On, optional observation O2n
Output: Observation of the reconstructed state o(u,.)
if then
F + CFE(o (o)) 01]
u1 Solver uo: F) return o(u1)
else
nn/2 + OP(o(u)), o, o( (unn22 + Reconstrruc(o((uu) On/2 On]
if O2n presen/ then
03n/2 OP O 02m)
on OP 03//2|
else
03n/2 absent
end
o( + Reconstrrccto((u..)) on, 03n/2|
return o(u,) end
The prediction refinement scheme, listed in Algorithm and illustrated in Figure 3c, solves this problem by re-evaluating existing predictions whenever the simulation progesses in time. Not all predictions need to be updated, though, and an update to prediction at a finer time scale can depend on a sequence other predictions. The prediction refinement algorithm that achieves this in an optimal form is listed in Appendix B. While the resulting execution order is difficult to follow for longer sequences with more than n = 8 frames, we give an overview of the algorithm by considering the prediction for time fn/2- After the first framm prediction o (tttn/)) of the n-frame sequence is made by OP,. the prediction refinement algorithm calls itself recursively until all frames up to frame n// are reconstructed from the CFE and the solver. The center prediction is then updated using OP,/2 for the next smaller time scale compared to the previous prediction. The call of OPn/2 also depends on (t(3nn/4). which was predicted using OP../- After half of the remaining distance to the center is reconstructed by the solver, the center prediction t,/2 is updated again, this time by the OP,, OP,/4. including all prediction dependencies. Hence, the center prediction continually refined every time the temporal distance between the latest reconstruction and the prediction halves, until

---

Published as a conference paper at ICLR 2020
(a) F(t) = 0
(b) CFE chain
(c) Supervised
(d) Diff. Phys.
(e) Ground Truth
1.0
to
0.5
t+
0.0
16
32 0
16 x (samples)
32 0
16 x (samples)
32 0
16 x (samples)
32 0
16 x (samples)
32
x (samples)
Figure 4: Trajectories for an example control task using Burger's equation. Initial and target states plotted with thick dashed lines in red and blue, respectively. Inferred states shown solid lines. (a) Natural evolution. (b) Reconstruction using a CFE chain. (c,d) Reconstructions using our hierarchical predictor-correcto scheme. (e) Ground-truth trajectory generated with constant force.
the reconstruction reaches that frame. This way, all final predictions (t,) are conditioned on the reconstruction of the previous state u(t;-1) and can therefore account for all previous deviations.
The prediction refinement scheme requires the same number of force inferences but an increased number of OP evaluations compared to the simpler algorithms. With total of 3n log, (n) OP evaluations (see Appendix B), it is of the same complexity, (()). In practice, this refinement scheme incurs only small overhead in terms of computation, which is outweighed by the significant gains in quality of the learned control function.
6 RESULTS
We evaluate the capabilities of our method to learn to control physical PDEs in three different test environments of increasing complexity. We first target simple but nonlinear 1D equation, for which we present an ablation study to quantify accuracy. We then study two-dimensional problems: an incompressible fluid and fluid with complex boundaries and indirect control. Full details are given in Appendix D. Supplemental material containing additional sequences for all of the tests can be downloaded from https:///ge in.tum...tt. ww..tum.eeeeee in..um.de/pubbicationss. httppppbccaaiinn/22-ii-i-------e
Burger's equation. Burger's equation is nonlinear PDE that describes the time evolution of single field, (LeVe,ue, 1992). Using Eq. it can be written as
Ou 8P
Ou 8ºu dz
P u, Ox dz2
Ox
Examples of the unperturbed evolution are shown in Figure 4a. We let the whole state be observable and controllable, i.e. (() = u(!), which implies that o* can always be reached exactly.
The results of our ablation study with this equation are shown in Table 1. The table compares the resulting forces applied by differently trained models when reconstructing a ground-truth sequence (Figure 4c). The variant denoted by CFE chain uses a neural network to infer the force without any intermediate predictions. With supervised loss, this method learns approximate single step well. However, for longer sequences, results quickly deviate from an ideal trajectory and diverge because the network never learned account for errors made previous steps (Figure 4b). Training the network with the objective loss (Eq. 4) using the differentiable solver greatly increases the quality of the reconstructions. On average, applies only 34% of the force used by the supervised model as it learns to correct the temporal evolution of the PDE model.
Next, we evaluate variants of our predictor-corrector approach, which hierarchically predicts in- termediate states. Here, the CFE is implemented as F(t.) = (o u(4)))Att. Unlike the simple CFE chain above, training with the supervised loss and staggered execution produces stable (albeit jittering) trajectories that successfully converge to the target state (Figure 4c). Surprisingly, this supervised method reaches almost the same accuracy as the differentiable CFE, despite not hav- ing access to physics- based gradients. However, employing the differentiable physies loss greatly

---

Published as a conference paper at ICLR 2020
Table 1: Quantitative reconstruction evaluation using Burger's equation, avg. for 100 examples.
Execution scheme
Training loss
Force f | dt
Inference time (ms)
CFE chain
Supervised
83.4 +: 2.0
0.024 0.013
CFE chain
Dift. Physics
28.8 + 0.8
4004
Staggered
Supervised
34.3 + 1.1
1.15
Staggered Refined
Diff. Physies
15.3
1.15 0.19
Diff. Physics
14.2 0.7
3.05 + 0.37
lterative optim. (60 iter.) Iterative optim. (300 iter.)
Diff. Physies
15.3
52.7 = 2.1
Diff. Physics
10.2 - 1.9
264.0 3.0
(a)
(b)
Figure 5: Example reconstructed trajectory from (a) the natural flow test set and (b) the shape test set. The target state o" shown on the right.
improves the reconstruction quality, producing solutions that are hard distinguish from ideal tra- jectories (Figure 4d). The prediction refinement scheme further improves the accuracy, but the differences to the staggered execution are relatively small as the predictions of the latter are already very accurate.
Table also lists the results of classic shooting-based optimization applied to this problem. To match the quality of the staggered execution scheme, the shooting method requires around 60 optimization steps. These steps are significantly more expensive to compute, despite the fast convergence. After around 300 iterations, the classic optimization reaches optimal value of 10.2 and the loss stops decreasing. Starting the iterative optimization with our method as an initial guess pushes the opti- mum slightly lower to 10.1. Thus, even this relatively simple problem shows the advantages of our learned approach.
Incompressible fluid flow. Next, we apply algorithm to two-dimensional fluid dynamics problems, which are challenging due to the complexities of the governing Navier-Stokes equa- tions (Batchelor, 1967). For a velocity field v, these can be written
P(u, Ve) - Ve V Vp,
subject the hard constraints V-v 0 and Vxp- = where p denotes pressure and the viscosity. In addition, we consider a passive density that moves with the fluid via 0p/80 0t = - Vp. We set be hidden and p to be observable, and allow forces to be applied to all of v.
We run our tests on 128 grid, resulting in more than 16,000 effective continuous control param- eters. We train the OP and CFE networks for two different tasks: reconstruction of natural fluid flows and controlled shape transitions. Example sequences are shown in Figure 5 and quantita- tive evaluation, averaged over 100 examples, given in Table 2. While all methods manage to approximate the target state well, there considerable differences in the amount of force applied. The supervised technique exerts significantly more force than the methods based on the differen- tiable solver, resulting in jittering reconstructions. The prediction refinement scheme produces the smoothest transitions, converging to about half the loss of the staggered, non-refined variant.
We compare method to classic shooting algorithms for this incompressible flow problem. While a direct shooting method fails to converge, a more advanced multi-scale shooting approach still re- quires 1500 iterations obtain level of accuracy that model achieves almost instantly. In
8

---

Published as a conference paper at ICLR 2020
Table 2:
comparison of methods in terms of final cost for (a) the natural flow setup and (b) the shape transitions. The initial distribution is sampled randomly and evolved to the target state.
Exccution
Loss
Force LF
Obs. L
b) Force LF b) Obs. L
Staggered
Supervised
243 1 11
1.53 0.23
n/a
n/a 0.331 0.134 0.126 = 0.010
Staggered
Diff. Physics
22.6 1.1
0.64 0.08
89
Refined
Diff. Physics
11.7 0.6
0.88 = 0.11
7 =
Figure 6: Indirect control sequence. Obstacles are marked white, control regions light blue. The white arrows indicate the velocity field. The domain enclosed in solid box with an open top.
addition, our model successfully learns solution manifold, while iterative optimization techniques essentially start from scratch every time. This global view leads our model to more intuitive solu- tions and decreases the likelihood of convergence to undesirable local minima. The solutions of our method can also be used as initial guesses for iterative solvers, as illustrated Appendix D.4. We find that the iterative optimizer with initial guess converges to solutions that require only 57.4% of the force achieved by the iterative optimizer with default initialization. This illustrates how the more global view of the learned solution manifold can improve the solutions of regular optimization runs.
Splitting the task into prediction and correction ensures that intermediate predicted states phys- ically plausible and allows us to generalize to new tasks. For example, we can infer transitions involving multiple shapes, despite training only individual shapes. This is demonstrated in Ap- pendix D.2.
Incompressible fluid with indirect control. The next experiment increases the complexity of the fluid control problem by adding obstacles to the simulated domain and limiting the area that can be controlled by the network. An example sequence in this setting is shown in Figure 6. As before, only the density p is observable. Here, the goal is to move the smoke from its initial position near the center into one the three "buckets" at the top. Control forces can only be applied in the peripheral regions, which are outside the visible smoke distribution. Only by synchronizing the 5000 continuous control parameters can a directed velocity field constructed the central region.
We first infer trajectories using trained CFE network and predictions that move the smoke into the desired bucket in straight line. This baseline manages to transfer 89% 2.6% of the smoke into the target bucket. Next we enable the hierarchical predictions and train the OPs. This version manages to maneuver 99. 99.22% + 0.15% of the smoke into the desired buckets while requiring 19.1% + 1.0% less foree.
For comparison, Table 3 also lists success rate and execution time for a direct optimization. Despite only obtaining low success of 82%, the shooting method requires several orders of magnitude longer than evaluating our trained model. Since all optimizations are independent of each other, some find better solutions than others, reflected in the higher standard deviation. The increased number of free parameters and complexity of the fluid dynamics to be controlled make this problem intractable for the shooting method, while our model can leverage the learned representation infer solution very quickly. Further details are given in Appendix D.3.
7 CONCLUSIONS
We have demonstrated that deep learning models in conjunction with a differentiable physics solver successfully predict the behavior of complex physical systems and learn to control them. The in-

---

Published as a conference paper at ICLR 2020
Table 3: Comparison of different methods the task of moving distribution of smoke into the target region by applying forces outside the region.
Method
Optimized quantity
Inside target (%)
Inference time (ms)
Straight trajectory
CFE
89.5 + 2.6
31.46 0.20
Staggered predictions
CFE, OP,
99.22 0.15
67.40 0.20
Iterative optim.
Control velocity
82.1
266.5 10
troduction of hierarchical predictor-correctt architecture allowed the model to learn to reconstruct long sequences by treating the physical behavior on different time scales separately.
We have shown that using differentiable solver greatly benefits the quality of solutions since the networks can learn how their decisions will affect the future. In our experiments, hierarchical infer- ence schemes outperform traditional sequential agents because they easily learn to plan ahead.
To model realistic environments, we have introduced observations to our pipeline which restrict the information available to the learning agent. While the PDE solver still requires full state information to run the simulation, this restriction does not apply when the agent is deployed.
While do not believe that learning approaches will replace iterative optimization, shows that it is possible to learn representations of solution manifolds for optimal control trajectories using data-driven approaches. Fast inference is vital in time-critical applications and can also be used in conjunction with classical solvers to speed up convergence and ultimately produce better solutions.
method
ACKNOWLEDCEMENTS
This work was supported in part by the ERC Starting Grant realFlow (BRC-2015-SGG-637014)
REFERENCES
Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. TensorFlow: system for large- scale machine learning. In Symposium on Operating Systems Design and Implementation, 2016.
Brandon Amos and Zico Kolter. OptNet: Differentiable optimization layer in neural networks. In ICML, 2017.
Yohai Bar-Sinai, Stephan Hoyer, Jason Hickey, and Michael P. Brenner. Learning data-driven dis- cretizations for partial differential equations. Proceedings of the National Academy of Sciences, 116(31), 2019.
George K. Batchelor. An Introduction to Fluid Dynamics. Cambridge University Press, 1967.
Peter W Battaglia, Jessica B Hamrick, and Joshua B Tenenbaum. Simulation an engine of physical scene understanding. Proceedings of the National Academy of Sciences, 110(45), 2013.
Thomas R. Bewley. Flow control: challenges for new renaissance. Progress in Aerospace Sciences, 37, 2001.
Hans Georg Bock and Karl-Josef Plitt. A multiple shooting algorithm for direct solution of optimal control problems. IFAC Proceedings Volumes, 17(2), 1984.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, and Skye Wanderman-Milne. JAX: Composable transformations of Python+NumPP programs. GitHub, 2018
John Bradshaw, Matt Kusner, Brooks Paige, Marwin HS Segler, and José Miguel Hernández- Lobato. A generative model for electron paths. In ICLR, 2019.
10

---

Published as a conference paper at ICLR 2020
Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differ- ential equations. In Advances in Neural Information Processing Systems, 2018.
Alexandre Joel Chorin. The numerical solution of the Navier-Stokes equations for an incompressible fluid. Bulletin of the American Mathematical Society, 73(6). 1967.
Richard Courant and David Hilbert. Methods of Mathematical Physics, volume II: Partial Differen- tial Equations. Interscience Publishers, 1962.
Filipe de Avila Belbute-Peres, Kevin Smith, Kelsey Allen, Josh Tenenbaum, and Zico Kolter. End-to-end differentiable physics for learning and control. In Advances in Neural Information Processing Systems, 2018.
Jonas Degrave, Michiel Hermans, Joni Dambre, and Francis wyffels. A differentiable physics engine for deep learning in robotics. Frontiers in Neurorobotics, 13, 2019.
Patrick E Farrell, David A Ham, Simon W Funke, and Marie E Rognes. Automated derivation of the adjoint of high-level transient finite element programs. SIAM Journal on Scientific Computing, 35(4), 2013.
Ciprian Foias, George R Sell, and Roger Temam. Inertial manifolds for nonlinear evolutionary equations. Journal of Differential Equations, 73(2), 1988.
David A Fournier, Hans J Skaug, Johnoel Ancheta, James Ianelli, Arni Magnusson, Mark N Maun- der, Anders Nielsen, and John Sibert. Ad model builder: Using automatic differentiation for statistical inference highly parameterized complex nonlinear models. Optimization Methods and Software, 27(2), 2012.
Michael B Giles and Niles A Pierce. An introduction to the adjoint approach to design. Flow, Turbulence and Combustion, 65, 2000.
Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl. Neural message passing for quantum chemistry. In ICML, 2017.
Gene H Golub and Charles F Van Loan. Matrix Computations. JHU press, 2012.
Andreas Griewank and Andrea Walther. Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. SIAM, 2nd edition, 2008
Andreas Griewank, David Juedes, and Jean Utke. Algorithm 755: ADOL-C: package for the automatic differentiation of algorithms written in C/C++. ACM Transactions on Mathematical Software (TOMS), 22(2), 1996.
Charles R Hargraves and Stephen W Paris. Direct trajectory optimization using nonlinear program- ming and collocation. Journal of Guidance, Control, and Dynamics, 10(4), 1987.
Francis Harlow and Eddie Welch. Numerical calculation time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, 8(12), 1965.
Carsten Hartmann, Juan C Latorre, Wei Zhang, and Grigorios A Pavliotis. Optimal control of multiscale systems using reduced-order models. Journal of Computational Dynamics, 1(2), 2014.
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recog- nition. CVPR, 2016.
Jun-Ting Hsieh, Shengjia Zhao, Stephan Eismann, Lucia Mirabella, and Stefano Ermon. Learning neural PDE solvers with convergence guarantees. In ICLR, 2019.
Yuanming Hu, Jiancheng Liu, Andrew Spielberg, Joshua B Tenenbaum, William Freeman, Jia- jun Wu, Daniela Rus, and Wojciech Matusik. ChainQueen: A real-time differentiable physical simulator for soft robotics. ICRA, 2019.
John Ingraham, Adam Riesselman, Chris Sander, and Debora Marks. Learning protein structure with differentiable simulator. In ICLR, 2019.

---

Published as a conference paper at ICLR 2020
Antony Jameson. Aerodynamic design via control theory. Journal of Scientific Computing, 3(3), 1988.
Dieter Kraft. On converting optimal control problems into nonlinear programming problems. In Computational Mathematical Programming. 1985.
Lubor Ladicky, SoHyeon Jeong, Barbara Solenthaler, Mare Pollefeys, and Markus Gross. Data- driven fluid simulations using regression forests. ACM Transactions Graphics (TOG), 34(6), 2015.
Horace Lamb. Hydrodynamics. Cambridge University Press, 1932.
Daniel B Leineweber, Irene Bauer, Hans Georg Bock, and Johannes P Schlöder. An efficient mul- tiple shooting based reduced SQP strategy for large-scale dynamic process optimization. Part I: Theoretical aspects. Computers Chemical Engineering, 27(2), 2003.
Randall J. Le Veque. Numerical Methods for Conservation Laws. Springer, 1992.
Yunzhu Li, Jiajun Wu, Russ Tedrake, Joshua B Tenenbaum, and Antonio Torralba. Learning particle dynamics for manipulating rigid bodies, deformable objects, and fluids. In ICLR, 2019.
Junbang Liang, Ming C. Lin, and Vladlen Koltun. Differentiable cloth simulation for inverse prob- lems. In Advances in Neural Information Processing Systems, 2019.
Julia Ling, Andrew Kurzawski, and Jeremy Templeton. Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. Journal of Fluid Mechanics, 807, 2016.
Jacques Louis Lions. Optimal Control of Systems Governed by Partial Differential Equations. Springer-Veriie Verlag, 1971.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. PDE-net: Learning PDEs from data, In ICML, 2018.
Dougal Maclaurin, David Duvenaud, and Ryan Adams. Autograd: Effortless gradicnts in Numpy. In ICML Workshops, 2015.
Antoine McNamara, Adrien Treuille, Zoran Popović, and Jos Stam. Fluid control using the adjoint method. ACM Transactions on Graphics (TOG), 23(3), 2004.
Arthur Mensch and Mathieu Blondel. Differentiable dynamic programming for structured prediction and attention. In ICML, 2018
Jeremy Morton, Antony Jameson, Mykel J Kochenderfer, and Freddie Witherden. Deep dynamical modeling and control of unsteady fluid flows. In Advances in Neural Information Processing Systems, 2018.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Fei-Fei Li, Josh Tenenbaum, and Daniel L. Yamins. Flexible neural representation for physics prediction. In Advances in Neural Information Processing Systems, 2018.
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-performan deep learning library. In Advances in Neural Information Processing Systems, 2019.
Olivier Pironneau. On optimum design in fluid mechanics. Journal of Fluid Mechanics, 64(1), 1974.
Lev Semenovich Pontryagin. Mathematical Theory of Optimal Processes. John Wiley, 1962.
Jovan Popovic, Steven M. Seitz, Michael A. Erdmann, Zoran Popovic, and Andrew P. Witkin. In- teractive manipulation of rigid body simulations. In SIGGRAPH, 2000.
William H Press, Saul Teukolsky, William T Vetterling, and Brian P Flannery. Numerical Recipes. Cambridge University Press, 3rd edition, 2007.
12

---

Published as a conference paper at ICLR 2020
Maziar Raissi, Alireza Yazdani, and George Karniadakis. Hidden fluid mechanics: A -Nviie--Stokee informed deep learning framework for assimilating flow visualization data. arXiv:1808.04333 2018.
Olaf Ronneberger, Philipp Fischer, and Thomas Brox. U-net:
Convolutional networks for biomedi- cal image segmentation. In Medical Image Computing and Compuler-Assisted Intervention, 2015.
Connor Schenck and Dieter Fox. SPNets: Differentiable fluid dynamics for deep neural networks.
In Conference on Robot Learning. 2018.
Gordon D Smith. Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press, 1985.
Jos Stam. Stable Fluids. In SIGGRAPH, 1999.
Edriss S Titi. On approximate inertial manifolds to the Navier-Stokes equations. Journal of Mathe- matical Analysis and Applications, 149(2), 1990.
Seiya Tokui, Ryosuke Okuta, Takuya Akiba, Yusuke Nitani,, Toru Ogawa, Shunta Saito, Shuji Suzuki, Kota Uenishi, Brian Vogel, and Hiroyuki Yamazaki Vincent. Chainer: A deep learn- ing framework for accelerating the research cycle. In Conference Knowledge Discovery and Data Mining (KDD), 2019.
Jonathan Tompson, Kristofer Schlachter, Pablo Sprechmann, and Ken Perlin. Accelerating Eulerian fluid simulation with convolutional networks. In ICML, 2017.
Mare Toussaint, Kelsey Allen, Kevin Smith, and Joshua B Tenenbaum. Differentiable physics and stable modes for tool-use and manipulation planning. In Robotics: Science and Systems, 2018.
Adrien Treuille, Antoine McNamara, Zoran Popovié, and Jos Stam. Keyframe control of smoke simulations. ACM Transactions on Graphics (TOG), 22(3), 2003.
Adrien Treuille, Andrew Lewis, and Zoran Popović. Model reduction for real-time fluids. ACM Transactions Graphics (TOG), 25(3), 2006.
Benjamin Ummenhofer, Lukas Prantl, Nils Thürey, and Vladlen Koltun, Lagrangian fluid simulation with continuous convolutions. In ICLR, 2020.
Bart van Merriënboer, Dan Moldovan, and Alexander Wiltschko. Tangent: Automatic differentiation using source-code transformation for dynamically typed array programming. In Advances in Neural Information Processing Systems, 2018.
Hermann von Helmholtz. Über integrale der hydrodynamischen gleichungen, welcher der wirbel- bewegungen entsprechen. Journal für Mathematik Bd. LV. Heft, 1:4, 1858.
Paul Werbos. Backwards differentiation in AD and neural nets: Past links and new opportunities. In Automatic Differentiation: Applications, Theory, and Implementations. Springer, 2006.
You Xie, Erik Franz, Mengyu Chu, and Nils Thuerey, tempoGAN: temporally coherent, volumet- ric GAN for super-resolution fluid flow. ACM Transactions Graphics (TOG), 37(4), 2018.
Kemin Zhou, John C. Doyle, and Keith Glover. Robust and Optimal Control. Prentice Hall, 1996.
13

---

Published as a conference paper at ICLR 2020
A IMPLEMENTATION DETAILS OF THE DIFFERENTIABLE PDE SOLVER
Our solver is publicly available at https://gitub..com/ com/tum com/tum-bss/Piilloo pbs/PhiFlow, licensed as MIT. It is imple- mented via existing machine learning frameworks to benefit from the built-in automatic differentia-- tion and to enable tight integration with neural networks.
For the experiments shown here used the popular machine learning framework Tensor- Flow (Abadi et al., 2016). However, our solver is written in a framework-independeer way and also supports PyTTorch (Paszke et al., 2019). Both frameworks allow for low-level NumPy-like implementation which is well suited for basic PDE building blocks. The following paragraphs out- line how we implemented these building blocks and how they can be put together to solve the PDEs shown in Section 6.
Staggered grids. Many of the experiments presented in Section use PDEs which track velocities. We adopt the marker and-cell marker-and-cell method (Harlow & Welch, 1965), storing densities in regular grid and velocity in a staggered grid. Unlike regular grids, where all components are sampled at the cen- ters of grid cells, staggered grids sample vector fields in staggered form. Each vector component sampled in the center of the cell face perpendicular to that direction. This sampling allows for an exact formulation of the divergence of staggered vector field, decreasing discretization errors in many cases. On the other hand, it complicates operations that combine vector fields with regular fields such as transport or density-dependent forces.
We use staggered grids for the velocities in all of our experiments. The buoyaney operation, which applies an upward force proportional to the smoke density, interpolates the density to the staggered grid. For the transport, also called advection, of regular or staggered fields, we interpolate the staggered field to grid cell centers or face centers, respectively. These interpolations are implemented in TensorFlow using basic tensor operations, similar to the differential operators. We implemented all differential operators that act on vector fields support staggered grids well.
Differential operators. For the experiments outlined in this paper, we have implemented the fol- lowing differential operators:
Gradient scalar fields in any number dimensions, V
Divergence of regular and staggered vector fields in any number of dimensions,
Curl of staggered vector fields in 2D and 3D, V
Laplace of scalar fields in any number of dimensions, V
All differential operators are local operations, i.e. they only act on small neighbourhood of grid points. In the context of machine learning. this suggests implementing them as convolution opera- tions with a fixed kerne.. Indeed, all differential operators can be expressed this way and we have implemented some low-dimensional versions of them using this method.
This method does, however, scale poorly with the dimensionality of the physical system the convolutional kernels pick up a large number of zeros, thus wasting computations. Therefore, we express n-dimensional differential operators using basic mathematical tensor operations.
Consider the gradient computation in ID, which results in staggered grid. Each resulting value is
(Vx); = 24 I;-1+
assuming the result is staggered at the lower faces each grid cell. This operation can imple- mented as ID convolution with kernel (-1, 1) as a vector operation which subtracts the array from itself, shifted by one element.
In low-dimensional setting, the convolution operation will be faster as it is highly optimized and can be executed on GPUs with one call. In higher dimensions, however, the vector-based version is faster and more practical because it avoids unnecessary computations and be coded in dimension-inen -nndeeednnnndeeeeeentt fashion. Both convolutions and basic mathematical operations are supported by all common machine learning frameworks, eliminating the need implement custom gradient functions.
Advection. PDEs containing material derivatives can be solved using an advection step Advect
14

---

Published as a conference paper at ICLR 2020
which moves each value of field f in the direction specified by vector field We implement the advection with semi- Lagianggann step (Stam, 1999) that looks back in time and supports regular and staggered vector fields.
To determine the advected value of grid cell or face Furget first v is interpolated to that point. Then the origin location is determined by following the vector backwards in time,
Torigin
Jiurget
Trarget
The final value is determined by linearly interpolating between the neighbouring grid cells around Forigin- All of these operations are, again, implemented using basic mathematical operations. Hence, gradients can be provided by the framework.
Poisson problems. Incompressible fluids, governed by the Navier-Stokes equations, are subject to the hard constraints V e = 0 and V x = O where p denotes the pressure. A numerical solver can achieve this by finding p such that these constraints are satisfied. This step is often called Chorin Projection, or Helmholtz decomposition, and is closely related to the fundamental theorem of vector calculus (von Helmholtz, 1858: Chorin, 1967). On grid, solving for is equal to solving a Poisson problem, i.e. system of N linear equations, Ap V where N is the total number of grid cells. The (N N) matrix A is sparse and its entries are located at predictable indices.
We numerically solve this Poisson problem with a conjugate gradient (CG) algorithm (Golub & Van Loan, 2012) that iteratively approximates p. Since hundreds of CG steps typically need to be performed for each Poisson solve, it is unfeasible to unroll this chain of iterations, and store all intermediate results in memory.
To ensure that the automatic differentiation chain is not broken, we instead solve the adjoint problem. For the pressure solve operation, the matrix A is symmetric and positive-definite This causes the adjoint problem to have the same mathematical form (McNamara et al., 2004) the original prob- lem. Therefore we implement the gradient for the pressure solve by performing a pressure solve on the gradient. We believe that this is good example of leveraging the methodology of adjoint method optimizations (Giles & Pierce, 2000; Treuille et al., 2006; Pontryagin, 1962) within deep learning. With this formalism we arrive at differentiable solver framework that closely integrates numerical methods for forward problems with support for inverse problems such as deep learning via the adjoint method.
Solving Burger's equation and the Navier-Stokee equations. Using the basic building blocks out- lined above, solving the PDEs becomes relatively simple. Burger's equation involves an advection and diffusion term which we evaluate independently.
Solver[uj [u]
Diffuse (Advect(u, ul]
where we explicitly compute Diffuse [u] + 2 with viscosity v. The advection is semi- Lagrangian with fack-tracing as described above.
Solving the Navier-Stoke equations, typically comprises of the following steps:
Transporting the density, p < Advect|p, v]
Transporting the velocity, v f- Advect|u, v]
Applying diffusion if the viscosity is v > O.
Applying buoyancy force, +- G p with buoyancy direction 0
Enforcing incompressibility by solving for the pressure, p <- Solve[Ap V v]. then
v n Vp
These steps are executed in this order to advance the simulation forward in time.
B COMPLEXITY OF EXECUTION SCHEMES
The staggered execution scheme recursively splits a sequence length n into smaller sequences, as depicted in Fig. 3b and Fig. 7a for 8. With each level of recursion depth, the sequence length
15

---

Published as a conference paper at ICLR 2020
(a) Staggered execution
(b) Prediction refinement
Fg
Time
to
Timet
Fg
(P(
OP(8)
(4(4)
OP(4)
(P(4)
(2)2)
(2(2)
OP(2)
Sal
Sal
OP
Sol
OP(2)
Sol
So
(P(4)
So
OP(4)
OP (2)
OP(2)
(P(2)
Sol
So
(P(2)
Sol
OP(2)
Sol
So
Sol
(P)
(2)
OP(2)
S0
(P(22)
Sol
OP(2)
Sol
So
Figure 7: OP and CFE+Solver (Sol) executions for sequence of length 8 performed by (a) the staggered execution scheme, (b) the prediction refinment scheme. The execution order is from top to bottom.
cut in half and twice as many predictions need to be performed. The maximum depth depends on the sequence length tn to and the time steps At performed by the solver,
f.
to
dmax
logg
- 1.
At
Therefore, the total number of predictions, equal to the number of OP evaluations, is
dras
Nop = 1 1
+ n/2
2-
2 = m-
0=
The prediction refinement scheme performs To understand the number of OP evaluations, we need to consider the recursive algorithm Reconstruct [uo, on: 02n]. listed in Alg that reconstructs sequence partial sequence of frames. For the first invocation, the last parameter O2n is absent, but for subsequences, that is not necessarily the case. Each invocation performs one OP evaluation O2n is absent, otherwise three. By counting the sequences for which this condition is fulfilled, we can compute the total number of network evaluations be
predictions,
be seen in Fig. 7b.
dies
Nop = 3
2 (o)g(n) = 3n 2 (n) -3
k==
C NETWORK ARCHITECTURES AND TRAINING
All neural networks used in this work are based on modified U-net architecture (Ronneberger et al., 2015). The U-net represents a typical multi-level convolutional network architecture with skip connections, which we modify by using residual blocks (He et al., 2016) instead of regular convolutions for each level. We slightly modify this basic layout for some experiments.
The network used for predicting observations for the fluid example is detailed Tab. 4. The input to the network are two feature maps containing the current state and the target state. Zero-padding is applied to the input, so that all strided convolutions do not require padding. Next, five residual blocks executed in order, each decreasing the resolution (1/2, 1/4, 1/8, 1/16, 1/32) while increasing the number of feature maps (4, 8, 16. 16, 16). Each block performs a convolution with kernel size and stride 2, followed by residual blocks with kernel size and symmetric padding. Inside each
16

---

Published as a conference paper at ICLR 2020
Table 4: Layers comprising the observation predictor network used in the direct fluid control exper- iment.
Layer
Resolution
Feature Maps
Input
128
2
Pad
159
2
Strided convolution + 2x Residual block 79
4
Strided convolution + 2x Residual block 39
8
Strided convolution 2x Residual block 19
16
Strided convolution + 2x Residual block
16
Strided convolution 2x Residual block
4
16
3x Residual block
16
Upsample Concatenate
32
Convolution + 2x Residual block
8
16
Upsample Concatenate
16
32
Convolution + 2x Residual block
16
16
Upsample Concatenate Convolution + 2x Residual block
32
24
32
16
Upsample Concatenate
64
20
Convolution 2x Residual block
64
16
Upsample Concatenate Convolution
128
18
128
1
block, the number of feature maps stays constant. Three more residual blocks are executed on the lowest resolution of the bowtie structure, after which the decoder part of the network commences, translating features into spatial content.
The decoder works as follows:
Starting with the lowest resolution, the feature maps are upsampled with linear interpolation. The upsampled maps and the output of the previous block of same reso- lution are then concatenated. Next, a convolution with 16 filters, kernel size of and symmetric padding, followed by two more residual blocks, is executed. When the original resolution is reached, only one feature map is produced instead of 16, forming the output of the network.
Depending on the dimensionality of the problem, either ID or 2D convolutions are used. The net- work used for the indirect control task is modified in the following ways: (i) It produces two output feature maps, representing the velocity (Vz, vy). (ii) Four feature maps of the lowest resolution (4x4) fed into dense layer producing four output feature maps. These and the other feature maps concatenated before moving to the upsampling stage. This modification ensures that the receptive field of the network is the whole domain.
networks were implemented TensorFlow (Abadi et al., 2016) and trained using the ADAM optimizer an Nvidia GTX 1080 Ti. We use batch sizes ranging from 4 to 16. Supervised training of all networks converges within a few minutes, for which we iteratively decrease the learning rate from 10-3 10-5 We stop supervised training after few epochs, comprising between 2000 and 10.000 iterations, as the networks usually converge within a fraction of the first epoch.
For training with the differentiable solver, start with decreased learning rate of 10-4 since the backpropagation through long chains is more challenging than training with a supervised loss. Optimization steps are also considerably more expensive since the whole chain needs be executed, which includes a forward and backward simulation pass. For the fluid examples, optimization step takes 1-2 seconds complete for the 2D fluid problems. We let the networks run about 100.000 iterations, which takes between one and two days for the shown examples.
D DETAILED DESCRIPTION AND ANALYSIS OF THE EXPERIMENTS
In the following paragraphs, give further details the experiments of Section 6.
17

---

Published as a conference paper at ICLR 2020
(a) = 0 (b) =c. (c) CFE Sup. (d) CFE Diff. (e) Stg. Sup. (f) Stg. Diff. (g) Refined
0
u
(a) F(()= (b) F(t) (c) CFE Sup. (d) CFE Diff. (e) Stg. Sup. (f) Stg. Diff.
(g) Refined
0
(a) F(t) 0 (b) F(t) (c) CFE Sup. (d) CFE Diff. (e) Stg. Sup. (f) Stg. Diff.
(g) Refined
0
-2
Figure 8: Trajectories for example control tasks using Burger's equation. The initial state is plot- ted in red, the target state in blue. (a) Natural evolution, with constant force, (c) CFE with supervised loss, (d) CFE with differentiable physics loss, (e) Re- construction with supervised loss, (f) Reconstruction with differentiable physics loss and staggered execution, (g) Reconstruction with differentiable physics loss and prediction refinement.
(b) Ground truth trajectory generated
D.1 BURGER'S ' EQUATION
For this experiment, we simulate Burger's equation (Eq. 6) on one-dimensional grid with 32 sam- ples over a course of 32 time steps. The typical behavior of Burger's equation in exhibits shock waves that move in + -r direction for u(x) > 0 or u(x) 0, respectively. When opposing waves clash, they both weaken until only the stronger wave survives and keeps moving. Examples shown in Figs. 4a and 8a.
All 32 samples are observable and controllable, i.c. o(t) = u(t). Thus, we can enforce that all trajectories reach the target state exactly by choosing the force for the last step to be
u((n-1))
F(tn-1)=
At
To measure the quality of a solution, it is therefore sufficient consider the applied force h" [F(1] which detailed for the tested methods in Table 1.
Network training. Both for the CFE chains as well as for the observation prediction models, we use the same network architecture, described in Appendix C. We train the networks on 3600 randomly generated scenes with constant driving forces, (t) = const. The examples are initialized with two Gaussian waves of random amplitude, size and position, set to clash in the center. In each time step, a constant Gaussian force with the same randomized parameters is applied to the system to steer it away from its natural evolution. Constant forces have larger impact the evolution than temporally varying forces since the effects of temporally varying forces can partly cancel out over time. The ground truth sequence therefore be regarded near-perfect but not necessarily
18

---

Published as a conference paper at ICLR 2020
optimal trajectory. Figs. 4d and 8b display such examples. The same trajectories, without any forces applied, are shown in sub-figures (a) for comparison.
We pretrain all networks (OPs or CFE, depending the method) with supervised observation loss.
L00 = OP (o)(4) ((;)) प. (4 11 2
(8)
The resulting trajectory after supervised training for the CFE chain shown in Figure 4b and Fig- ure Sc. For the observation prediction models, the trajectories are shown in Figure 4c and Figure Sc.
After pretraining, we train all OP networks end-to-end with our objective loss function (see Eq. 4), making use of the differentiable solver. For this experiment, we choose the mean squared difference for the observation loss function:
- (o(((t)))
0²12. -012²
We test both the staggered execution scheme and the prediction refinement scheme, shown in Fig- ure 8f and Figure 8g.
Results. Table 1 compares the resulting forces inferred by different methods. The results are aver- aged over a set of 100 examples from the test set which is sampled from the same distribution as the training set. The CFE chains both fail to converge to o". While the differentiable physies version manages to produce a u that resembles o", the supervised version completely deviates from an optimal trajectory. This shows that learning to infer the control force F(t:) only from u(t:), o" and t is very difficult as the model needs to learn to anticipate the physical behavior over any length of time.
Compared to the CFE chains, the hierarchical models require much less force and learn to converge towards o². Still, the supervised training applies much more force the system than required, the reasons for which become obvious when inspecting Figure 4b and Fig. 8c. While each state seems close the ground truth individually, the control oscillates undesirably, requiring counter-actions later in time.
The methods using the differentiable solver significantly outperform their supervised counterparts and exhibit an excellent performance that is very close the ground truth solutions in terms of required forces. On many examples, they even reach the target state with less force than was applied by the ground truth simulation. This would not be possible with the supervised loss alone, but by having access to the gradient-based feedback from the differentiable solver, they can learn to find more efficient trajectories with respect to the objective loss. This allows the networks to learn applying forces in different locations that make the system approach the target state with less force.
Figure 4e and Fig. show examples of this. The ground truth applies the same force in each step, thereby continuously increasing the first sample u(x = 0), and the supervised method tries to imitate this behavior. The governing equation then slowly propagates u(r 0) in positive x direction since u(x = 0) > O. The learning methods that use a differentiable solver make use of this fact by applying much more force F(x- 0) 0 at this point than the ground truth, even overshooting the target state. Later, when this value had time to propagate to the right, the model corrects this overshoot by applying negative force F(r = 0) 0. Using this trick, these models reach the target state with up to 13% less force than the ground truth on the sequence shown in Figure 4.
Figure 9 analyzes the variance of inferred forces. The supervised methods often fail to properly converge to the target state, resulting in large forces in the last step, visible as second peak in the supervised CFE chain. The formulation of the loss (Eq. 3) suppresses force spikes. In the solutions inferred by our method, the likelihood of large forces falls off multi-exponentially as a consequence. This means that large forces are exponentially rare, which the expected behavior given the 1.2 regularizer from Eq. 3.
We also compare our results to a single-shooting baseline which is able to find near-optimal solutions at the cost of higher computation times. The classic optimization uses the ADAM optimizer with learning rate of 0.01 and converges after around 300 iterations. To reach the quality of the staggered prediction scheme, it requires only around 60 iterations. This quick convergence can be explained by the relatively simple setup that is dominated by linear effects. Therefore, the gradients are stable, even when propagated through many frames.
19

---

Published as a conference paper at ICLR 2020
10a
CFE chain (Supervised)
Staggered execution (Supervised)
10
CFE chain (Diff. Phys.)
Prediction Refinement (Diff. Phys.)
10a
Counts
10a
101
10a
0.0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6
Force per point
Figure 9: Histogram comparing the frequency of force strengths of different methods, summed over 100 examples from the Burger's experiment.
The computation times, shown Tab. 1, were recorded on a single GTX 1080 Ti. We run 100 examples in parallel to reduce the relative overhead caused by GPU instruction queuing. For the network-based methods, we average the inference time over 100 runs. We perform 10 runs for the optimization methods.
D.2 INCOMPRESSIBLE FLUID FLOW
The incompressible Navier-Stokes equations model dynamics of fluids such as water or air, which develop highly complex and chaotic behavior. The phenomenon of turbulence is generally seen as one of the few remaining fundamental and unsolved problems of classical physics. The challenging nature of the equations indicates that typically a very significant computational effort and a large number of degrees of freedom are required to numerically compute solutions. Here, we target an incompressible two-dimensional gas with viscosity v, described by the Navier-Stokes equations for the velocity field v. We assume a constant fluid density throughout the simulation, setting Py const. The gas velocity controllable and, according Eq. 1, we set
Vv) = -(v V)u Vp
PS
subject to the hard constraints V v = 0 and V x = 0. For our experiments, we target fluids with low viscosities, such as air, and v apply numerical diffusion that is on average higher than the targeted one. For fluids with a larger viscosity, the Poisson solver outlined above computing p could be used to implicitly solve vector-valued diffusion equation for v.
the equation above as the transport steps implicitly
However, incorporating significant amount of viscosity would make the control problem easier to solve for most cases, as viscosity suppresses small scale structures in the motion. Hence, in order to create challenging environment for training our networks, we have but minimal amount of diffusion in the physical model.
In addition to the velocity field v, we consider a smoke density distribution which moves passively with the fluid. The evolution of is described by the equation 0p/Ot - V.. We treat the velocity field as hidden from observation, letting only the smoke density be observed, i.e. o(t) = p(t). We (v,p) to write the system
stack the two fields
PDE, compatible with Eq. 1.
20

---

Published as a conference paper at ICLR 2020
(a) Supervised, staggered execution
oP, 6
u
(b) Diff. Physics, staggered execution
OP,
u
(c) Diff. Physics, prediction refinement
oP.
r
Figure 10: Reconstruction of an example natural flow sequence from the test set. The predictions o" are plotted above the reconstructions u. The target state o is shown in the last column of the predictions.
21

---

Published as a conference paper at ICLR 2020
For the OP and CFE networks, we use the 2D network architecture described in Appendix C. Instead of directly generating the velocity update in the CFE network for this problem setup, we make use of stream functions (Lamb, 1932). Hence, the CFE network outputs vector potential Φ of which the curl V x + is used as a velocity update. This setup numerically simplifies the incompressibility condition of the Navier-Stokes equations but retains the same number of effective control parame- ters.
Datasets. We generate training and test datasets for two distinct tasks: flow reconstruction and shape transition. Both datasets have a resolution of 128 x 128 with the velocity fields being sam- pled in staggered form (see Appendix A). This results in over 16.000 effective continuous control parameters that make up the control force F(t:) for each step i.
The flow reconstruction dataset is comprised of ground-truth sequences where the initial states (Pu, vo) are randomly sampled and then simulated for 64 time steps. The resulting smoke density is then taken to be the target state, o" = p" prim (f6a). Since we use fully convolutional networks for both CFE and OPs, the open domain boundary must be handled carefully. If smoke was lost from the simulation, because it crossed the outer boundary, neural network would see the smoke simply vanish unless it was explicitly given the domain size input. To avoid these problems, we run the simulation backwards in time and remove all smoke from Po that left the simulation domain.
For the shape transition dataset, we sample initial and target states po and p. by randomly choosing a shape from a library containing ten basic geometric shapes and placing it at a random location inside the domain. These can then be used for reconstructing sequences of any length n. For the results on shape transition presented in section 6, we choose n = 16 because all interesting behavior can be seen within that time frame. Due the linear interpolation used in the advection step (see Appendix A), both and v smear out over time. This numerical limitation makes it impossible to match target states exactly in this task as the density will become blurry over time. While we could generate ground-truth sequences using classical optimizer, we refrain from doing so because (i) these trajectories are not guaranteed be optimal and (ii) we want to see how well the model can learn from scratch, without initialization.
Training. We pretrain the CFE on the natural flow dataset with a supervised loss,
L CFE (t) = (uu(t) + F(t) (t)
where (L) denotes the velocity from ground truth sequences. This supervised training alone con- stitutes good loss for the CFE as it only needs to consider single-step intervals At while the OPs handle longer sequences. Nevertheless, we found that using the differentiable solver with an obser- vation loss
I B, (0') - B, (Solver(u + CFE(u, o']]) 12
further improves the accuracy of the inferred force without sacrificing the ground truth match. Here 1
B, (x) denotes a blur function with kernel of the form 1+2/7. The blur helps make the gradients
smoother and creates non-zero gradients in places where prediction and target do not overlap. During training, we start with a large radius of r = 16 Ar for B, and successively decrease it to r = 2 Ar. We choose a such that LF and Lgg of the same magnitude when the force loss spikes (see Fig. 15).
After the CFE is trained, we successively train the OPs starting with the smallest time scale. For the OPs, we train different models for natural flow reconstruction and shape transition, both based on the same CFE model. We pre-train all OPs independently with supervised observation loss before jointly training them end-to-end with objective loss function (Eq. 4) and the differentiable solver find the optimal trajectory. We use the OPs trained with the staggered execution scheme as initialization for the prediction refinement scheme. The complexity of solving the Navier-Stokes equations over many time steps this example requires such a fully supervised initialization step. Without it. this setting is so non-linear that the learning process does not converge to good so- lution. Hence, it illustrates the importance of combining supervised and unsupervised (requiring differentiable physics) training for challenging learning objectives.
A comparison of the different losses is shown in Fig. 10. The predictions, shown in the top rows of each subfigure, illustrate the differences between the three methods. The supervised predictions, especially the long-term predictions (central images), are blurry because the network learns to aver- age over all ground truth sequences that match the given initial and target state. The differentiable
22

---

Published as a conference paper at ICLR 2020
(a) Predicted trajectory
(b) Reconstructed trajectory
Reconstruction of multiple shapes using prediction refinement. The CFE is trained the natural flow dataset and OPs are trained on single-shape transitions. The predictions of all shapes added and passed the CFE one prediction.
Figure 11:
physics solver largely resolves this issue. The predictions are much sharper but the long-term pre- dictions still do account for short-term deviations. This can be seen in the central prediction of Fig. 10b which shows hints of the target state o", despite the fact that the actual reconstruction u cannot reach that state at that time. The refined prediction, shown in subfigure (c), is closer to u since it is conditioned on the previous reconstructed state.
In the training data, we let the network transform one shape into another at a random location. The differentiable solver and the long-term intuition provided by our execution scheme make it possible to train networks that can infer accurate sequences of control forces. In most cases, the target shapes closely matched. As networks infer sequences over time, we refer readers to the supplemental material (https://ge.in wum.de/publications//000 (itpiiooos/222-0--ooooll)) which contains animations of additional sequences.
Generalization to multiple shapes. Splitting the reconstruction task into prediction and correction has the additional benefit of having full access the intermediate predictions o². These model real states of the system so classical processing or filter operations can be applied to them as well. We demonstrate this by generalizing our method m shapes that evolve within the same domain. Figure 11 shows an example of two weakly-interacting shape transitions. We implement this by executing the OPs independently for each transition E 21, while inferring the control force Ft) on the joint system. This is achieved by adding the predictions of the smoke density before passing
the CFE network, O' Σi=1 of The resulting force is then applied to all sequences individually so that smoke from one transition does not end up in another target state. Using this scheme, we can define start and end positions for arbitrarily many shapes and let them evolve together.
Evaluation of foree strengths The average force strengths are detailed in Tab. 2 while Figure 12 gives more detailed analysis of the force strengths. As expected from using a L2 regularizer on the force, large values are exponentially rare in the solutions inferred from our test set. None of the hierarchical execution schemes exhibit large outliers. The prediction refinement requires the least amount of force to match the target, slightly ahead of the staggered execution trained with the same loss. The supervised training produces trajectories with reduced continuity that result in larger forces being applied.
D.3 INCOMPRESSIBLE FLUID WITH INDIRECT CONTROL
As fourth test environment, we target case with increased complexity, where the network does not have the means anymore to directly control the full fluid volume. Instead, the network can only apply forces the peripheral regions, with a total of more than 5000 control parameters per step. The obstacles prevent fluid from passing through them and the domain is enclosed with solid boundaries from the left, right and bottom. This leads to additional hard constraints and interplays between constraints in the physical model, and such provides interesting and challenging test case for our method. The domain has three target regions (buckets) separated by walls at the top of the domain, into which volume of smoke should be transported from any position in the center
23

---

Published as a conference paper at ICLR 2020
Staggered execution (Supervised)
Staggered execution (Diff. Phys.)
10a
Prediction Refinement (Diff. Phys.)
10a
Counts
10a
101
0
1
2
4
5
6
Force per point
Figure 12: Histogram comparing the frequency of force strengths applied in the direct fluid control experiment the natural flow dataset, summed over 100 examples.
Figure 13: Two reconstructed trajectories from the test set of the indirect smoke control problem.
24

---

Published as a conference paper at ICLR 2020
0.8
0.4
50
150
250
350
450
550
Iterative optimization of the indirect incompressible fluid control problem. The graph shows the fraction of smoke that ends up in the correct bucket vs number of optimization steps, averaged over 10 examples.
Figure 14:
part. Both initial position and the target bucket are randomized for our training set of 3600 examples and test set of 100 examples. Each sequence consists of 16 time steps.
In this case the control indirect since the smoke density lies outside the controlled area at all times. Only the incompressibility condition allows the network to influence the velocity outside the controlled area. This forces the model consider the global context and synchronize a large number of parameters to create desired flow field. The requirement of complex synchronized force fields makes generating reliable training data difficult, as manual or random sampling is unlikely to produce directed velocity field in the center. We therefore skip the pretraining process and directly train the CFE using the differentiable solver, while the OP networks are trained as before with 2 Az.
To evaluate how well the learning method performs, we measure how much the smoke density ends up inside the buckets and how much force applied in total. For reference, replace the observation predictions with an algorithm that moves the smoke towards the bucket in a straight line. Averaged over 100 examples from the test set, the resulting model manages put 89% + 2.6% of the smoke into the target bucket. In contrast, the model trained with our full algorithm moves 99.22% + 0.15% of the smoke into the target buckets while requiring 19.1% + 1.0% less force.
We also compare our method to an iterative optimization which directly optimizes the control ve- locities. We use the ADAM optimizer with a learning rate of 0.1. Despite the highly non-linear setup, the gradients are stable enough to quickly let the smoke flow in the right direction. Fig. 14 shows how the trajectories improve during optimization. After around 60 optimization steps, the smoke distribution starts reaching the target bucket in some examples. Over the next 600 iterations, it converges to a a configuration in which 82.1 + 7.3 of the smoke ends up in the correct bucket.
D.4 COMPARISON TO SHOOTING METHODS
We compare the sequences inferred by our trained models to classical shooting optimizations using our differentiable physics solver to directly optimize F't) with the objective loss L (Eq. 4) for a single input. We make use of stream functions (Lamb, 1932), as in the second experiment, to ensure the incompressibility condition is fulfilled. For this comparison, the velocities of all steps are initialized with normal distribution with μ = 0 and o 0.01 so that the initial trajectory does not significantly alter the initial state, u(()= & u(to).
We first show how a simple single-shooting algorithm (Zhou et al., 1996) fares with our Navier- Stokes setup. When solving the resulting optimization problem using single-shooting, strong ar- tifacts in the reconstructions can be observed, as shown in Figure 17a. This undesirable behavior stems from the nonlinearity of the Navier-Stokes equations, which causes the gradients Au > to become noisy and unreliable when they recurrently backpropagated through many time steps. Unsurprisingly, the single-shooting optimizer converges to a undesirable local minimum.
As single-shooting is well known to have problems with non-trivial problem settings, employ multi-scale shooting (MS) method (Hartmann et al., 2014). This solver first computes the trajectory on coarsely discretized version of the problem before iteratively refining the discretization. For
25

---

Published as a conference paper at ICLR 2020
(a) Observation loss
(b) Force loss
4.000e+4
240
4.000e+4
200
160
2.000e14
120
1.000c14
80.0
40.0
0.000
2.000k
4.000k
6.000k
0,000
2,000k
4,000k
8.000k
Figure 15: Mean convergence curves of the adjoint method optimization for 100 shape transitions.
(a) Observation loss
(b) Force loss
24.0
200
20.0
160
16.0
120
12.0
80.0
0.00
2.000k
4.000K
0.000
2.000k
4.C00k
6.000k
Figure 16:
Mean convergence curves of the adjoint method optimization for 100 shape transitions, taking the reconstruction from the refinement scheme as initial guess.
the first resolution, we use 1/16 of the original width and height which both reduces the number of control parameters and reduces nonlinear effects from the physics model. By employing an exponential learning rate decay, this multi-scale optimization converges reliably for all examples. We use the ADAM optimizer to compute the control variable updates from the gradients of the differentiable Navier-Stokes solver.
An averaged set of representative convergence curves for this setup is shown in Figure 15. The objective loss (Eq. 4) shown its decomposed state as the sum of the observation loss I. shown in Figure 15a, and the force loss LF, shown in Figure 15b. Due to the initialization of all velocities with small values, the force loss starts out small. For the first 1000 iteration steps, 4G dominates which causes the system to move towards the target state o". This trajectory is not ideal, however, as more force than necessary is applied. Once observation loss and force loss are of the same magnitude, the optimization refines the trajectory to use less force.
We found that the trajectories predicted by our neural network based method correspond to perform- ing about 1500 steps with the MS optimization while requiring less tuning. Reconstructions of the same example are compared in Figure 17. Performing the MS optimization up to this point took 131 seconds on GTX 1080 Ti graphies card for a single 16-frame sequence while the network inference ran for 0.5 seconds. For longer sequences, this gap grows further because the network inference time scales with (n)) This could only be matched if the number of iterations for the MS optimization scaled with O(1), which is not the case for most problems. These tests indicate that our model has successfully internalized the behavior of a large class of physical behavior, and can exert the right amount of force to reach the intended goal. The large number of iterations required for the single-case shooting optimization highlights complexity of the individual solutions.
Interestingly, the network also benefits from the much more difficult task to learn whole manifold of solutions: comparing solutions with similar observation loss for the MS algorithm and our net- work, the former often finds solutions that are unintuitive and contain noticeable detours, e.g., not taking a straight path for the density matching examples of Fig. 5. In such situations, our network benefits from having to represent the solution manifold, instead of aiming for single task optimiza- tions. As the solutions are changing relatively smoothly, the complex task effectively regularizes the inference of new solutions and gives the network more global view. Instead, the shooting optimiza-
26

---

Published as a conference paper at ICLR 2020
(a) Classical optimization
(b) Classical optimization with multi-resolution
(c) Diff. physics, prediction refinement
Figure 17: Example reconstruction of shape transition. (a) Direct shooting optimization, 2300 iterations, (b) multi-scale shooting optimization, 1500 iterations, (c) output of our neural network based method with prediction refinement. Our model infers the shown solution in single pass, and generalizes to large class of inputs.
tions have to purely rely on local gradients for single-shooting or manually crafted multi-resolution schemes for MS.
Our method can also be employed to support the MS optimization by initializing it with the velocities inferred by the networks. In this case, shown in Figure 16, both L; and Lp decrease right from the beginning, similar to the behavior in Figure 15 from iteration 1500 on. The reconstructed trajectory from the neturalletworr-baaa described above is not necessary.
method is close to the optimum that the mustutresolution approach
D.5 ADDITIONAL RESULTS
In Fig. 18, we provide visual overview of a sub-set of the sequences that can be found in the supplemental materials. contains 16 randomly selected reconstructions for each of the natural flow, the shape transitions, and the indirect control examples. In addition, the supplemental mate- rial, available at https://ge.in.tum wecppubbicatiion//220-ccee unsupervised, staggered, and refined versions of our approach.
highlights the differences between
27

---

Published as a conference paper at ICLR 2020
(a) Natural flow
(b) Shape transition
Figure 18: Five additional sequences from the test sets of the natural flow and shape transition setups. The first nine frames contain frames from our reconstruction. The far right image shows the target.
28


Reviews:
[{"review_id": "rkeA9bdX5r", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "[Summary]\n\nThis paper proposes to combine deep learning and a differentiable PDE solver for understanding and controlling complex nonlinear physical systems over a long time horizon. The method introduces a predictor-corrector scheme, which employs a hierarchical structure that temporally divides the problem into more manageable subproblems, and uses models specialized in different time scales to solve the subproblems recursively.\n\nFor dividing the problem into subproblems, they use an observation predictor network to predict the optimal center point between two states. To scale the scheme to sequences of arbitrary length, the number of models scales with O(log N). For each subproblem, the authors propose to use a corrector network to estimate the control force to follow the planned trajectory as close as possible. \n\nThey have compared their method with several baselines and demonstrated that the proposed approach is both more effective and efficient in several challenging PDEs, including the incompressible Navier-Stokes equations.\n\n[Major Comments]\n\nPredicting the middle point between two states for modeling the dynamics via deep neural networks is not new, but I did not know any other works that use this idea for controlling PDEs.\n\nI like the idea of splitting the control problem into a prediction and a correction phase, which leverages the power of deep neural networks and also incorporates our understanding of physics. The introduction of the hierarchical structure alleviates the problem of accumulating error in single-step forwarding models and significantly improves the efficiency of the proposed method. The videos for fluid control in the supplement materials also convincingly demonstrate the effectiveness of the technique.\n\nI still have a few questions regarding the applicability and the presentation of the paper. Please see the following detailed comments.\n\n[Detailed Comments]\n\nIn Section 3, the authors claim that their model \"is conditioned only on these observables\" and \"does not have access to the full state.\" However, the model requires a differentiable PDE solver to provide the gradient of how interactions affect the outcome. These seem to contradict each other. Doesn't the solver require full-state information to predict the behavior of the system?\n\nRelated to the previous question, how can we make use of the differentiable PDE solver if we are uncertain or unknown of the underlying physics, i.e., partially observable scenarios.\n\nThe algorithm described in Section 5 seems to be the core contribution of this work. Instead of describing the algorithm in words, I think it would make it more clear if the authors can add an algorithm block in the main paper. It would also be better if the authors can include a few sentences describing the algorithm in the abstract to inform the readers of what to expect.\n\nFigure 4 is a bit confusing, and it would be better if the authors can include the label for the x-axis. Besides, in the caption, the authors said that they show \"the target state in blue.\" However, there are a lot of blue lines in the figure, and it is hard to know, at first glance, which one of them is the target.\n\nIn Table 1, the bottom two methods are using the same execution scheme and training loss, but the results are different. Is there a typo? Also, it would be better to bold the number that has the best performance.", "extra_scores": {"experience_assessment": "I have published one or two papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "rJgsfMWpYH", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "In this paper, the authors outline a method for system control utilizing an \"agent\" formed by two neural networks and utilizing a differentiable grid-based PDE solver (assuming the PDE describing the system is known). The agent is split into a control force estimator (CFE) which applies a force to advance the state of the controlled system, and an observation predictor (OP) which predicts the trajectory needed to reach the target state. The objective is to reach the target state with minimal total amount of applied force. The order of CFE and OP calls is discussed and the importance of keeping the trajectory predictions conditioned on the actual previous state of the system, so that errors from previous steps can be taken into account.\n\nThree application examples are discussed: Burger's equation (1D), and incompressible flow (2D) with direct and indirect control. In all cases the proposed scheme of \"prediction refinement\" leads to better or comparable results than standard iterative optimization, and is much more computationally efficient (at inference time, not taking into account the cost of training).\n\nThe paper presents an interesting mix of neural networks and traditional PDE solvers for system control, and I vote for acceptance. An additional advantage of the paper is the authors' promise to open source their differentiable PDE solver implemented in TensorFlow, which should make it easy for others to build upon their work. The text is easy to read, but quite verbose, with many of the technical details relegated to the (sizeable) appendices. I would recommend trying to trim it down where possible (for instance, the description of the U-nets could be more compact, perhaps in table form; 2nd paragraph of the background section seems a bit out of context and could probably be omitted, etc).\n\nQuestions and suggestions for improvements:\n\n* What form of L_o^* and alpha was used in all the experiments?\n* It looks like pretraining was used for all cases except for the most challenging one with indirect control. Was it truly necessary for the simpler experiments? \n* Improve naming consistency. It looks like \"differentiable physics\" and \"differentiable solver\" are used for the same thing in different parts of the paper. My recommendation would be to use the latter term everywhere.\n* How many time steps are used in the indirect control experiment?\n* IIUC, the optimization of Eq. 3 is always done end-to-end. Have any experiments been done to estimate how many time steps can be reliably handled by the proposed procedure before the optimization problem becomes too hard?\n", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}, {"review_id": "Sye8KBAe9B", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": null, "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "## Summary\n\nThe authors propose a method for training physical systems whose behavior is governed by partial differential equations. They consider situations where only partial observations available, and where control is indirect.\n\nIndirectly controlling physical systems is a very important problem with applications throughout engineering. In fact, much of the field of robotics can be described in these terms. The authors employ a number of interesting methods, including a predictor-corrector framework and the adjoint sensitivity method for differentiating through differential equation solvers.\n\nThe paper is generally very clear, organized and well written. There are only a few places where I think clarification is needed (see detailed comments below). I also have a few questions about the losses and training procedure. On the whole, I think the paper is inventive, well-written and potentially very impactful. I think it would be a great addition to ICLR.\n\n\n## Clarifications\n\n* Page 4: I found the statement \"an agent trained in supervised fashion will then learn to average over the modes instead of picking one of them\" a little confusing. Could you clarify the reasoning here?\n* Page 5: I think the description of predictor-corrector could be clearer. In particular, I found the phrase \"the correction uses o(t + ∆t) to obtain o(t + ∆t)\" unclear.\n* Page 8: Could you add a description of what is observable to the body of the paper (I see it is included in the supplement)?\n* Page 8: I think ∇p needs to be divided by density in your NS equation, right?\n* Page 8 - 9: Is there a limit on the size of the force that can be applied at any point? I know the total force is penalized, but what about the maximum force applied at any point?\n\n\n## Losses and Training\n\n* I think \"differentiable physics\" losses need a more detailed explanation in the body of the paper.\n* In the supplement, it is defined using B_r, but I don't think B_r is defined.\n* It seems like the differential physics loss requires a differential solver (in this case, for Burger/Navier-Stokes). If I have understood this correctly, I think this needs to be discussed in the body of the paper. In particular, it would be nice to discuss what happens when the physics is a black box (i.e. we can interact with the system by applying control and observing, but we don't know the rules governing the physical system). Is this exactly when we are restricted to the \"supervised\" loss? Is there some middle ground? What if we had black box access to the exact physics, along with an approximate differentiable solver? This seems like a realistic scenario for e.g. large fluid flow scenarios.", "extra_scores": {"experience_assessment": "I have read many papers in this area.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper at least twice and used my best judgement in assessing the paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a method to control dynamical systems described by a partial differential equations (PDE). The method uses a hierarchical predictor-corrector scheme that divides the problem into smaller and simpler temporal subproblems. They illustrate the performance of their method on 1D Burger’s PDE and 2D incompressible flow.
The reviewers are all positive about this paper and find it well-written and potentially impactful. Hence, I recommend acceptance of this paper.
Decision: Accept (Spotlight)