{'fraction_non_alphanumeric': 0.03798848325599598, 'fraction_numerical': 0.012886453849164198, 'mean_word_length': 4.503846153846154, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatio-Temporal Stochastic Optimization: Theory and Applications to Optimal Control and Co-Design

There is a rising interest in Spatio-temporal systems described by Partial Differential Equations (PDEs) among the control community. Not only are these systems challenging to control, but the sizing and placement of their actuation is an NP-hard problem on its own. Recent methods either discretize the space before optimziation, or apply tools from linear systems theory under restrictive linearity assumptions. In this work we consider control and actuator placement as a coupled optimization problem, and derive an optimization algorithm on Hilbert spaces for nonlinear PDEs with an additive spatio-temporal description of white noise. We study first and second order systems and in doing so, extend several results to the case of second order PDEs. The described approach is based on variational optimization, and performs joint RL-type optimization of the feedback control law and the actuator design over episodes. We demonstrate the efficacy of the proposed approach with several simulated experiments on a variety of SPDEs.


Introduction
In many complex natural processes, a variable such as temperature or displacement has values that are time varying on a spatial continuum over which the system is defined. These spatio-temporal processes are typically described by Partial Differential Equations (PDEs) and are increasingly prevalent throughout the robotics community. Swarm robotics can be described by Reaction-advection-diffusion PDEs [1]. Robot navigation in crowded environments can be described by Nagumo-like PDEs [2]. Soft robotic limbs can be modelled as damped Euler-Bernoulli systems [3].
Some of the major control-related challenges of these systems include dramatic under-actuation, high system dimensionality, and the design and/or placement of distributed actuators over a continuum of potential locations. These systems often have significant time delay from a control signal, and can have several bifurcations and multi-modal instabilities. In addition, realistic representations of these systems are stochastic. Put together, control of spatio-temporal systems represents many of the largest current-day challenges facing the robotics and automatic control community. This paper addresses stochastic optimal control and co-design of Stochastic Partial Differential Equations (SPDEs) through the lens of stochastic optimization. We propose a joint actuator placement and policy network optimization strategy via episodic reinforcement that leverages inherent spatio-temporal stochasticity in the dynamics for optimization. The resulting stochastic gradient descent approach bootstraps off the widespread success of SGD methods such as ADAM for training Artificial Neural Networks (ANNs).
Among the main goals of this line of research is to establish capabilities for the eventual design and manufacture of soft-body robots. The behavior of such systems follows second order SPDEs. As such, while the proposed method is general to first and second order systems, we focus our mathematical formulation on second order SPDE.
Our Contributions. In this work we tackle the coupled challenge of policy optimization and actuator co-design for SPDEs. Our approach is founded on a general principle coming from thermodynamics that also has had success in stochastic optimal control literature [4] Free Energy ≤ Work − Temperature × Entropy. (1) We leverage this principle in order to derive a measure-theoretic loss function that utilizes exponential averaging over importance sampled system trajectories in order to choose network and actuator design parameters that simultaneously minimize state cost and control effort. This work builds off related work [5] by including actuator co-design and establishing needed theory for applicability to second order SPDEs. Specifically, we contribute the following: (i) A derivation of a joint policy optimization and actuator co-design architecture in Hilbert spaces for SPDEs (ii) A version of the Girsanov Theorem suitable for second-order SPDEs in a natural form (iii) A practical set of tools to extend related work to actuator co-design and to handle second-order SPDEs

Related Work
Control of SPDEs. There is a growing body of work that seek control of PDEs by immediately reducing them to a set of ODEs [6,7,8,9,10]. They do not consider stochasticity and typically use standard tools from finite-dimensional machine learning and control theory. Other approaches apply SPDE control theory, yet either treat linear systems or lack numerical results [11,12]. We build the proposed framework on [5], wherein the authors create a semi-model-free reinforcement learning framework for policy-based control of SPDEs. This related work does not consider actuator co-design and does not give a mathematical treatment suitable to second-order SPDEs.
Actuator Co-Design for PDEs. In [13], theoretical results for the optimal actuator location for minimum norm control of the stochastic heat equation are obtained, however, without a numerical demonstration. Optimal co-design techniques for actuator and sensor locations with H ∞ and H 2 design objectives were demonstrated on flexible structures [14,15,16], and the linearized Ginzburg-Landau equation [17,18], but are only applicable to linear systems. Similarly, symmetry properties and linear control theory are applied to chaotic systems in [19]. Optimal actuator and sensor placement for scalar control of the linear advection equation are demonstrated in [20,21] using Gramian based approaches which also requires system linearity. Similarly, [22] provides a level set method that promises scalability, yet relies on linear system Grammians. Conditions for the existence of optimal actuator placement for semilinear PDEs are obtained in [23]. Optimal actuator and sensor placement for stabilization of the nonlinear Kuramoto-Sivashinsky equation is demonstrated in [24], however they impose strong simplifying assumptions which limit their dimensionality.

Problem Formulation
This work proposes actuator co-design optimization for a large class of stochastic spatio-temporal systems represented as SPDEs. We describe these systems as evolving on time separable Hilbert spaces, where they are represented by infinite dimensional vectors and acted on by operators. We address a class of SPDEs that are of semi-linear form. Let H denote a separable Hilbert space with σ-field B(H) and probability space (Ω, F, P) with filtration F t , t ∈ [0, T ]. Consider the general semi-linear form of a controlled SPDE on H given by where X(t) ∈ H is the state of the system which evolves on the Hilbert space H, the linear and nonlinear measurable operators A : H → H and F (t, X) : R × H → H (resp.) are uncontrolled drift terms, Φ(t, X, x; Θ (k) ) : R × H × D → H is the nonlinear control policy parameterized by Θ (k) at the k th iteration, where D ⊂ R 3 is the domain of the finite spatial region, dW (t) : R → H is a spatio-temporal noise process, and G(t, X) is a nonlinearity that affects both the noise and the control. It is used to incorporate the effects of actuation on the field. The Hilbert spaces formulation given in eq. (2) is general in that any semi-linear SPDE can be described in this form by appropriately choosing the A and F operators. In this form, the spatio-temoral noise process dW (t) is a Hilbert space-valued Wiener process, which is a generalization of the Wiener process in finite dimensions. We include a formal definition of a Wiener process in Hilbert spaces for clarity [25, Section 4.1.1] be a complete orthonormal system for the Hilbert Space H. Let Q denote the covariance operator of the Wiener process W (t). Note that Q satisfies Qe i = λ i e i , where λ i is the eigenvalue of Q that corresponds to eigenvector e i . Then, W (t) ∈ H has the following expansion: where β j (t) are real valued Brownian motions that are mutually independent on (Ω, F , P).
This expansion in eq. (3) reveals how the Wiener process acts spatially. There are various forms of the Wiener process with different properties. We refer the interested reader to [25] for a more complete introduction. The proposed approach is derived for a special case of Wiener process called the Cylindrical Wiener process, defined as follows.
Definition 3.2. A Wiener process W (t) on H is called a Cylindrical Wiener process if the covariance operator Q is the identity operator I.
Note that for the Cylindrical Wiener process, the sum in eq. (3) is unbounded in H since λ j = 1, ∀j = 1, 2, . . . . This makes the Cylindrical Wiener process a challenging Wiener process to handle since it acts spatially everywhere with equal magnitude, in contrast to Wiener processes with covariance operators that are of trace class (i.e. wherein the expansion eq. (3) is finite). This type of spatio-temporal noise requires the assumption that the operators A , F (t, X), and G(t, X) satisfy properly formulated conditions given in [25,Hypothesis 7.2] to guarantee the existence and uniqueness of the F t -adapted weak solution X(t), t ≥ 0.
The nonlinear policy Φ(t, X, x; Θ (k) ) is a potentially time-varying policy that has explicit state dependence. Nonlinear, explicit state dependence allows for a feedback policy that can extract pertinent information from the state for control. In this work the nonlinear policy utilizes an ANN. Embedded in this function is also a spatial dependence. The dependence on x describes how the actuator is placed in the spatial domain. This approach encompasses cases where terms that parametrize how the actuators are shaped or sized are included in the nonlinear policy.
Many complex spatio-temporal systems are given by partial differential equations of second order in time. One such system is the simply supported stochastic Euler-Bernoulli equation with Kelvin-Voigt and viscous damping, which can describe the motion of a soft robotic limb. Formally, this is given by where the spatial region is one-dimensional, y(x, t) = y : R n ×R → R represents the vertical displacement of the beam, and all functional dependencies of the nonlinear policy Φ have been dropped since it has a different form in the PDE perspective. With the change of variables v := ∂ t y, this system has the typical second order matrix form where A 0 = ∂ xxxx without boundary conditions. Now, we lift this PDE into infinite dimensional Hilbert spaces. Define Y ∈ H as the Hilbert space analog of y(x, t), V ∈ H as the Hilbert space analog of v(x, t), and a variable Z on the direct product Hilbert space H 2 := H × H. Note that Z is a Hilbert space analog of a variable z(x, t) = [y(x, t) v(x, t)] ∈ R 2 . In Hilbert spaces, A 0 becomes an operator acting on H and 1 gets replaced by the identity operator I acting on H. Rewriting eq. (5) in Hilbert space semi-linear form yields where A : , G : H → H 2 is an operator representing how control and spatio-temporal noise enter the system G = [0; I], and dW (t) is a Cylindrical Wiener process on H. Note that the Hilbert space variables Y , V , and Z no longer have spatial dependence as the Hilbert space vectors capture the spatial continuum over which the problem is defined.
In order to derive a measure theoretic view of variational optimization, we require a change of measures coming from an appropriate Girsanov theorem for second order Hilbert space valued systems. We present the second order version of the Girsanov theorem here as an extension of existing formulations, where all norms and inner products are taken with respect to the Hilbert space H.
Let Ω be a sample space with a σ-algebra F. Consider the following H 2 -valued nonlinear stochastic processes where Z(0) =Z(0) = z 0 and W ∈ H is a Cylindrical Wiener process with respect to measure P. Moreover, let Γ be a set of continuous-time, infinite-dimensional trajectories in the time interval [0, T ]. Now the probability law of Z will be defined as L(Γ) := P(ω ∈ Ω|Z(·, ω) ∈ Γ) . Similarly, the law ofZ is defined asL(Γ) := P(ω ∈ Ω|Z(·, ω) ∈ Γ). Theñ where we have defined and assumed E P e Under the above assumption,Ŵ is a Cylindrical Wiener process with respect to a measure Q defined by The proof for this intermediate result can be found in [26,Theorem 10.14]. Now, using eq. (11), eq. (7) is rewritten as Notice that the SPDE in eq. (14) has the same form as eq. (8). Therefore, under the introduced measure Q and noise profileŴ , Z(·, ω) becomes equivalent toZ(·, ω). Conversely, under measure P, eq. (13) (or eq. (14)) behaves as the original system in eq. (7). In other words, eq. (7) and eq. (14) describe the same system on (Ω, F, P). From the uniqueness of solutions and the aforementioned reasoning, one has The result follows from eq. (12).
It follows that the Randon-Nikodym (RN) derivative between measures L(·) andL(·) of the different dynamical systems defined in eqs. (7) and (8) is given by Applying this to both the general semi-linear system eq. (2) and the second order semi-linear form of the Euler-Bernoulli system eq. (6) yields ψ := √ ρΦ(t, Z, x; Θ (k) ) and a Radon-Nikodym derivative as Throughout the rest of the manuscript, we refer to the terms in the Radon-Nikodym derivative many times. For convenience, they will be assigned functions N (Θ,

Spatio-Temporal Stochastic Optimization
The proposed measure theoretic framework is based on eq. (1) in the following form [27,28] where J = J(X) is an arbitrary state cost function. Relating eq. (17) to eq. (1), the metaphorical work and entropy describe a metaphorical energy landscape for which there is a minimizing measure. Sampling from this measure would simultaneously minimize state cost and the KL-divergence term, which is interpreted as control effort. The measure that optimizes eq. (17) is the so-called Gibbs measure It is not known how to sample directly from the Gibbs measure in eq. (18). Instead, variational optimization methods seek to iteratively minimize the controlled distribution's distance to the Gibbs measure [28,29,30]. Define the control policy and actuator co-design problem as Throughout experiments, the authors found that a joint optimization problem dramatically outperforms the split problem in eq. (19). To make this clear, define a new variableΘ := [Θ; x] , and with it the new joint variational optimization asΘ * = argmin Θ Expanding the KL divergence and performing importance sampling yieldŝ which is equivalent to minimizingΘ * = argmin Θ log dL dL dL * .
Performing importance sampling again yieldŝ The proposed iterative approach performs episodic reinforcement with respect to a loss function in order to optimize eq. (23). Define the loss function as Plugging eq. (18) and eq. (15) into eq. (24) yields whereJ ) is defined as and J(Z T 0 ) is a state cost evaluated over the state trajectory Z T 0 . For reaching tasks, J(Z T 0 ) is typically a weighted quadratic penalization of the 1-norm distance to the goal state.
This loss function compares sampled trajectories by evaluating them on the exponentiatedJ performance metric. The importance sampling terms N and P, which appear inJ add a quadratic control penalization term and a mixed control noise term. In the Loss function, they serve as weights for the exponentiated cost trajectories. For convenience, we denote the exponentiated cost term as Recall, that the nonlinear policy Φ is a functional mapping into Hilbert space H. This is kept general for derivation purposes, however it implies that the nonlinear policy affects controls each element of an infinite vector Z ∈ H. In a realistic context, there may not be control over how the actuation affects the field state Z. A more realistic, but less general representation refines the policy as where m(x) : D N → R N × H represents the effect of the actuation from N actuators on the infinitedimensional field. This function can take many forms, however in most applications it is either a Gaussian-like exponential with mean centered at the actuator locations, or an indicator function for actuation.
In eq. (27), ϕ(X; Θ (k) ) : H → R N is a policy network with N control outputs representing N distributed (or boundary) actuators. Note that as desired, the tensor contraction given on the right hand side of eq. (27) produces a vector in H. Splitting the actuation function from the control signal is also desired because we ultimately wish to use a finite input, finite output policy network for the function ϕ(X; Θ (k) ). The inner product terms become where M(x) = m(x)m(x) .

Approximate Discrete Optimization
Performing spatio-temporal stochastic optimization in Hilbert spaces as described above maintains generality of the resulting loss function. In addition, it allows us to make use of the spatio-temporal noise process. Finally, we can apply the spatial inner product integration scheme described in [5]. Despite this, computations on a digital computer require spatial discretization. As a result, the portion of this optimization procedure dedicated to actuator co-design becomes a discrete optimization problem. To see this more clearly, consider the 1-dimensional (1D) spatial continuum D = [0, 1] discretized into a 10 point 1D grid. Lets assume that an actuator is chosen to be placed at x = 0.25. Even though the actuation function m(x) may be Gaussian-like function, the majority of the actuation will be felt in between two grid points, namely 0.2 and 0.3. This problem is even more severe if the actuation function m(x) is the indicator function, as there will be no actuation exerted on the field irrespective of the control signal magnitude. Denote the number of spatial discretization points as J and the discretized problem domain grid asD composed of J 3 elements. The optimization problem becomes This formulation is an accurate representation, yet limits gradient flow from the loss function back to the actuator design parameters. In order to maintain these gradients, we approximate eq. (30) as follows. Define a one-to-one map S :D → Z + , where Z + denotes positive integers. Applying the forward and inverse mapping produces a gradient-based parameter update of the form where R(·) simply rounds to the nearest integer, γ Θ is the learn rate for the ANN parameters, γ x is the learn rate for the actuator design parameters, ∇ Θ denotes the gradient with respect to Θ, and ∇ x deontes the gradient with respect to x. This approach allows us to leverage well-known backprop-based algorithms such as ADA-Grad [31] and ADAM [32].

Algorithm and Network Architecture
As discussed previously, implementation of the above framework requires spatial and temporal discretization of the SPDEs discussed in the first section. With this in mind, we choose an ANN for our nonlinear for r = 1 to R do 5: for t = 1 to T do 6: dW t ← SampleN oise() 7: Z t ← P ropagate(Z t−1 , Θ (k) , dW t ) via eq. (6) 8: u t ← SparseF orwardP ass(Θ (k) , Z t ) Compute ∇ Θ (L) via backprop 17: Compute ∇ x (L) via backprop 18: (31) 19: x (k) ← GradientStep(L, γ x ) via eq. (32) 20: end for policy ϕ(Z; Θ (k) ). In this work we exclusively use Feed-forward Neural Networks (FNNs) for all of our experiments, and use physics-based models of each SPDEs to generate training data. Given that the proposed framework is semi-model-free, real system data can seamlessly replace the physics-based model as described in [5]. We only need prior knowledge of the flavor of noise and the actuator dynamics (i.e. m(x)). Also described is a sparse method for spatial integration that we apply here to maintain computational and memory efficiency.
In addition to obtaining gradients from any variant of Stochastic Gradient Descent (SGD) via backprop paths, GradientStep for actuator parameters also adds the current gradient to the previous gradient if the actuator location has not moved by the current gradient step. This is useful because when the actuators are near their optimal value, the gradients become very small, preventing the actuator locations from reaching optimal values. This could also be achieved by heuristically changing the learning rate γ x . Also, a quadratic cost term was added to L in order to penalize actuators from leaving the spatial domain, with a large coefficient to ensure this condition was rarely, if ever violated.
Note that there are separate learning rates for the policy and the actuators. This is because in practice the authors found that the optimization landscape is typically much more shallow for the actuator design than for the policy parameters. For most of the experiments, the actuator placement learn rate γ x was about 30 times larger than the policy network learn rate γ Θ .

Simulation Results and Discussion
We applied our approach to four simulated SPDE experiments to simultaneously place actuators and optimize a policy network. Each experiment used less than 32 GB RAM, and was run on a desktop computer with a Intel R Xeon(R) 12-core CPU with a NVIDIA GeForce GTX 980 GPU. Our code was written to operate inside a Tensorflow graph [34] to leverage rapid static graph computation, as well as sparse linear algebra operations used by SparseF orwardP ass. The first two experiments involved a reaching task, where the SPDEs is initialized at a zero initial condition over the spatial region, and must reach certain values at pre-specified regions of the spatial domain. The last two experiments involved a suppression task, where some non-zero initial condition must be suppressed on desired regions.
The data that was used for training was generated by a spatial central difference, semi-implicit time discretized version of each SPDE. These schemes are described in detail in [33,Chapter 3 & 10]. Each experiment had all actuators initialized by sampling from a random distribution on [0.4a, 0.6a], where a denotes the spatial size. For 3500 iterations of our algorithm, run times for the most complicated system-the Euler-Bernoulli equation-were about 15 hours.
The first experiment was a temperature reaching task on the 1D Heat equation with homogeneous, T = 0 Dirichlet boundary conditions, and is depicted in fig. 1. The task was to raise the temperature at regions specified in green to specified values depicted in the figure. The algorithm was run for 3000 iterations.
The next experiment was a velocity reaching task on the Burgers equation with non-homogenous Dirichlet boundary conditions, and is depicted in fig. 2. The Burgers equation has a nonlinear advection term, which produces an apparent rightward motion. The algorithm was run for 3500 iterations, and was able to take advantage of the advection for actuator placement in order to solve the task with lower control effort.
The heat equation is a pure diffusion SPDE, while the Burgers equation shares the diffusion term with the Heat equation with an added advection term. The results of the Heat and Burgers experiments show actuator locations that take advantage of the natural behavior of each SPDEs. In the case of the Heat equation, actuators are nearby the desired regions such that the temperature profile can reach a flat peak of the diffusion at the desired profile. In the case of the burgers equation, the advection pushes towards the right end of the space, thus forming a wave front that develops at the right end, but leaves the left end dominated by the diffusion term. This is again reflected in the placement of actuators. The first actuator is nearby the desired region just as the actuators in the Heat SPDE, while two of the actuators between the center and the right region are located to be able to control the amplitude and shape of the developing wave front so as to produce a flat peak that aligns with the desired region at the desired velocity. The central desired region is flanked on both sides by actuators that are nearly equidistant, in order to produce another desired flat velocity region at this location.
The third experiment was a voltage suppression task on the Nagumo equation with homogeneous Neumann boundary conditions, and is depicted in fig. 3. The task was to suppress an initial voltage on the left end, that without intervention propagates toward the right end, as shown by the uncontrolled trajectories. The Nagumo equation is composed of a diffusive term and a 3rd-order nonlinearity, making this equation the most challenging from a nonlinear control perspective. Despite this, our approach was able to simultaneously place actuators and provide control such that the task was solved. The algorithm was run for 2000 iterations, and demonstrates actuator placement optimization that takes advantage of the natural system behavior. This task was also the most challenging due to the significantly longer planning horizon of 3.5 seconds, as compared to the 1.0 second planning horizon of all the other experiments.
In order to validate our proposed approach, we compared the actuator locations that the algorithm found after optimization to the actuator locations that were hand placed by a human expert for the simulated experiments conducted in related work [5]. To have a valid comparison, we ran the IDVRL algorithm for  both sets of actuator locations. Figure 4 reports these results. The left figure shows that the state costs for each is almost identical. Note that the scale here is 10 −4 . The center figure shows the control signals for each actuator, for each method, and demonstrates that for almost identical state cost values, the control effort for each actuator with our approach is lower on average. The calculated average control signal magnitudes for IDVRL are 3.3 times higher than our method. The third plot shows the voltage profile at the final time. We hypothesize that the lower control effort is due to the control over the shape of the spatially propagating signal, enabling it to have a smoother transition into the desired region. While the penalty of this actuator placement is a higher variance on the desired region, the choice appears correct given the result.
The final task was an oscillation suppression task on the Euler-Bernoulli equation with Kelvin-Voigt damping given in eq. (4), and is depicted in fig. 5. As shown, the initial condition prescribes spatial oscillations, that then oscillate temporally. The second-order nature of the system creates offset and opposite oscillations in the velocity profile, that in turn produce offset and opposite oscillations in the   position profile. Without interference, the oscillations proceed over the entire time window. As shown on the right, our approach successfully suppresses these oscillations, which die out quickly under the given control policy. In this experiment, the actuators did not leave the initialized actuator placement region [0.4, 0.6] perscribed for all experiments.
The Euler-Bernoulli oscillation suppression task is in fact very challenging and complex. Producing a control signal at an actuator location that is in phase with the velocity oscillations will amplify the oscillations, leading to a divergence. The actuator location and control signal from the policy network must work in concert to produce a control signal out of phase with the velocity that matches its frequency, which is time varying due to control, as shown on the right side of fig. 5.
Each of the above experiments has its challenges and in most cases the spatio-temporal problem space produces a joint policy optimization and actuator co-design problem that is littered with local minima. These experiments demonstrate that the proposed approach can jointly optimize a policy network and actuator design. These results and the overall performance of the algorithm are indicative that this approach may enable actuator design on problem spaces where a human has little to no prior knowledge to rely on in attempting to solve the problem by hand.

Conclusion and Future Directions
This work presents a framework for joint policy optimization and actuator co-design. We contribute a novel mathematical tool required for a measure theoretic treatment of second order SPDEs, namely the change of measure and associated proof. We also provide a mechanism for handling discrete parameter optimization with a gradient-based approach. We demonstrate the resulting algorithm on four different SPDEs, each with their own challenges and complexities. The last of which is the Euler-Bernoulli SPDE which connects back with our goal of establishing capabilities for the further development of soft-body robotics.
The presented approach is a new way of performing optimization and can lead to many applications in soft robotics, soft materials, morphing and continuum mechanism. The results are encouraging to the authors. We plan to scale this approach to systems with higher spatial dimensionality (i.e. 2D and 3D), as well as investigate more complex forms of the above systems that are relevant to a future of soft robotics.

S1 Additional Information on Simulations
The following are details on each of our simulated experiments which will help in reproducing our results.
In every simulated experiment we used spatial central difference discretization and semi-implicit time discretization. Our time discretization size was ∆t = 0.01 throughout all experiments. All experiments initialized the actuator locations by sampling a uniform distribution on [0.4a, 0.6a] at the policy network initialization phase. All experiments utilized FNN for the nonlinear policy ϕ(h; Θ), with two hidden layers and each neuron applying ReLU activations. The number of neurons per layer matched the state dimension after discretization, which ended up being 64 for all experiments. All policy network weights were initialized with the Xavier initialization [35]. In every experiment the function m(x) was modeled as a Gaussian-like exponential function with the means co-located with the actuator locations. Each experiment added a penalization term if actuators left the spatial region, with a coefficient of 10 3 . Each experiment considered a cost function of the form where S is the spatial subregion on which the desired profile is defined. These desired spatial regions were different for each experiment and are detailed below.

S1.1 Heat SPDE Reaching Task
The controlled 1D heat SPDE in Hilbert space form with homogeneous Dirichlet boundary conditions is given by where is the thermal diffusivity parameter and was set to 1 for this experiment. The task is to achieve a desired temperature profile at three regions along the spatial domain of length a = 1.0 with three actuators. The spatial region was discretized into 64 equally spaced grid points. This experiment used a noise term ρ = 10. Here, the variance of the N = 3 actuator Gaussian-like functions was set to σ 2 = (0.1a) 2 for all actuators. In this experiment G(t, h) is the identity operator. The desired regions were defined by 78a, 0.82a]} is the right spatial region.
(S4) and the desired temperature was set to T = 1.0 for x ∈ S 1 ∪ S 3 , and T = 0.5 for x ∈ S 2 . The network was trained using the ADAM optimizer for 3000 iterations with 200 trajectory rollouts sampled from the Heat SPDE model per iteration. Each trajectory was 1.0 seconds long. The policy learn rate γ Θ was set to 10 −3 and the actuator location learn rate γ x was set to 3 × 10 −2 .

S1.2 Burgers SPDE Reaching Task
The controlled 1D Burgers SPDE in Hilbert space form with non-homogeneous Dirichlet boundary conditions is given by where the parameter is the viscosity of the medium. Equation (S5) considers a simple model of a 1D flow of a fluid in a medium with non-zero flow velocities at the two boundaries. The goal again here is to reach certain flow velocities at three specified spatial regions in the spatial domain of length a = 1.0. The spatial region was discretized into 64 equally spaced grid points. Here we consider G(t, h) as an identity operator. This experiment used a noise term ρ = 10. In this experiment the N = 5 Gaussian-like functions used a variance of σ 2 = (0.1a) 2 for every actuator. The scaling parameter κ was set to 100. The desired regions were given by eq. (S4). The desired values were given by h desired (t, x) = 2.0 m/s for x ∈ S 1 ∪ S 3 , which is at the sides, and h desired (t, x) = 1.0 m/s for x ∈ S 2 , which is in the central region.
The network was trained using the ADAM optimizer for 3500 iterations with 100 trajectory rollouts sampled from the Burgers SPDE model per iteration. Each trajectory was 1.0 seconds long. The policy learn rate γ Θ was set to 10 −3 and the actuator location learn rate γ x was set to 3 × 10 −2 .

S1.3 1D Nagumo SPDE Suppression Task
The 1D Nagumo equation in Hilbert space form with homogeneous Neumann boundary conditions is given by where the parameter α = −0.5 determines the speed of a wave traveling down the length a = 5.0 of the axon and = 1.0 determines the rate of diffusion. We consider G(t, h) as an identity operator. The goal of this task was to suppress the activation voltage that traveled rightward without control intervention. The N = 3 actuator Gaussian-like functions used a variance of σ 2 = (0.1a) 2 for each actuator. Here, the noise term was set to ρ = 10.0 The spatial domain was discretized using a grid of 64 points. The cost function coefficient κ was chosen as 10 −3 . The experiment used one desired spatial region S = [0.7a, 0.99a], which is at the right end. The network was trained using the ADAM optimizer for 2000 iterations with 100 trajectory rollouts sampled from the Nagumo SPDE model per iteration. Each trajectory was 3.5 seconds long. The policy learn rate γ Θ was set to 10 −3 and the actuator location learn rate γ x was set to 5 × 10 −2 .

S1.4 Euler-Bernoulli SPDE Suppression Task
The controlled 1D Euler-Bernoulli SPDE with homogeneous Dirichlet boundary conditions is given in the paper, but repeated here for completess ∂ tt y + ∂ xx ∂ xx y + C d ∂ xxt y + µ∂ t y = Φ + 1 √ ρ ∂ t W (t) y(t, 0) = y(t, a) = 0 y(0, x) = y 0 ∂ t y(0, x) = v 0 ∂ xx (t, 0) + C d ∂ xxt y(t, 0) = 0 ∂ xx (t, a) + C d ∂ xxt y(t, a) = 0 (S7) This is a second order system in time. Similar to a typical treatment of second order Stochastic Ordinary Differential Equation (ODE)s (SDEs), we define a velocity state v := ∂ t y. After rewriting (S7) and lifting it into Hilbert spaces, we have dZ = AZdt + G Φ(t, Z, x; Θ (k) )dt + 1 √ ρ dW (t) , y(t, 0) = y(t, a) = 0 where Z contains the both the displacement y and displacement velocity v, A : H 2 → H 2 is the linear operator A = [0 I; −A 0 − C d A 0 − µI], G : H → H 2 , and dW (t) is a Cylindrical Wiener process. The Kelvin-Voigt damping term C d was set to 10 −4 , and the viscous damping term µ was set to 10 −3 . This was to almost completely remove damping of the initialized oscillations throughout the entire time window. The initial displacement was set by the function As shown in the uncontrolled contour plot ( fig. 5 in the main paper), this initial condition generates oscillations over the spatial length a = 1.0 that evolve temporally. The position and velocity spaces were each spatially discretized into 32 equal points, for a total state space dimensionality of 64. In this experiment, the N = 8 actuator Gaussian-like functions used a variance of σ 2 = (0.2a) 2 . The noise term ρ was set to 1.0 for this experiment. The cost function coefficient κ was set to 3 × 10 −4 , and the desired spatial subregion was defined over the entire space S = [0, 1.0]. The network was trained using the ADAM optimizer for 3500 iterations with 250 trajectory rollouts sampled from the Euler-Bernoulli SPDE model per iteration. Each trajectory was 1.0 seconds long. The policy learn rate γ Θ was set to 10 −3 and the actuator location learn rate γ x was set to 5 × 10 −2 .