Optimal Prediction using Learning and Shape Optimization

Estimation is of fundamental importance in many disciplines. For example, it is important in climate science to predict the temperature changes in the Earth atmosphere (Park et al., 2019). Estimation is also used to find sea surface temperature in order to predict weather (Nielsen-Englyst et al., 2018). In economics, estimation is used to predict investment risks. In control theory, estimation is used to feedback from the state of the system at locations not accessible by sensors. Most commonly, estimation is used to reconstruct the state of a system modeled by partial differential equations such as the heat and wave equations.

In recent years, machine learning methods have been used to extract predictive models from large sets of data provided by sensors. These methods can further find a symbolic representation of the model (symbolic regression) and also discover the hidden physics. Among the machine learning methods, deep neural networks have gained considerable attention (LeCun et al., 2015). Evolutionary algorithms (Vaddireddy et al., 2020), compressive sensing (Candès & Wakin, 2008), and sparse optimization (Candès et al., 2008) are also machine learning algorithms used for logistic regression. In (Long et al., 2018), a PDE-net approach is proposed to predict the dynamics of systems modeled by partial differential equations (PDEs) from sensor data. The network is also able to uncover the underlying PDE model. In (Lu et al., 2018), it is shown that deep neural networks such as ResNet, PolyNet, FractalNet, and RevNet can be interpreted as different numerical discretizations of differential equations. This interpretation is further used to design new deep neural networks. In (Liu et al., 2010, 2013), training data for image restoration is used to create a learning-based PDE model for computer vision tasks. The problem is stated as an optimal control problem where the inputs and outputs are training images. The authors show the effectiveness of their model by numerical experiments on image denoising and inpainting problems in image restoration.

The shape optimization of sensors and actuators in the context of PDEs has been discussed in relatively few works. In (Privat et al., 2013), the optimal shape and position of an actuator for the wave equation in one spatial dimension are discussed. An actuator is placed on a subset $\omega\in[0,\pi]$ with a constant Lebesgue measure $L\pi$ for some $L\in(0,1)$. The optimal actuator minimizes the norm of a Hilbert Uniqueness Method (HUM)-based control; such control steers the system from a given initial state to zero state in finite time. In (Privat et al., 2017), the optimal actuator shape and position for linear parabolic systems are discussed. This paper adopts the same approach as in (Privat et al., 2013) but with initial conditions that have randomized Fourier coefficients. The cost is defined as the average of the norm of HUM-based controls. In (Kalise et al., 2018), optimal actuator design for linear diffusion equations has been discussed. A quadratic cost function is considered, and shape and topological derivatives of this function are derived. Numerical results show significant improvement in the cost and performance of the control. In (Edalatzadeh et al., 2019), the optimal shape of actuators for vibration control of flexible beams is investigated using a linear-quadratic performance index and shape optimization.

Optimal sensor design problems are in many ways similar to the optimal actuator design problem. In (Privat et al., 2014), optimal sensor shape design has been studied where the observability is being maximized over all admissible sensor shapes. Optimal actuator design problems for nonlinear distributed parameter systems have also been studied. In (Edalatzadeh & Morris, 2019), it is shown that under certain conditions on the nonlinearity and the cost function, an optimal input and actuator design exist, and optimality conditions are derived. Results are applied to the nonlinear railway track model as well as to the semi-linear wave model in two spatial dimensions. Numerical techniques to calculate the optimal actuator shape design are mostly limited to linear quadratic regulator problems, see, \eg, (Allaire et al., 2010). For controllability-based approaches, numerical schemes have been studied in (Münch & Periago, 2011; Münch, 2007; Münch & Periago, 2013).

In this paper, we combine shape optimization with learning-based PDE models to find the optimal shape of sensors, as well as an estimator of the full state of the PDE. In \crefsection:estimator_design, the neural-network predictor is formulated, which recovers the full PDE state from limited sensor data and predicts the state at the subsequent time step. \Crefsection:shape_optimization discusses the shape optimization of sensors. Simulations results are presented in \crefsection:simulation_results. Concluding remarks are given in \crefsection:conclusion.

Consider a first-order dynamical system with state $z(t)$, initial condition $z_{0}$, state space $X$, and linear or nonlinear operator $F$ on $X$:

In our applications, $F$ is a second-order partial differential (PDE) operator acting on functions in one or two spatial dimensions. The evolution of the state $z(t)$ in discrete time with sufficiently small time increment $\Delta t$ can be described, using for simplicity an explicit Euler method, as

with initial value $z_{0}$. For numerical computation, we also need to discretize in space. For simplicity of notation, we currently assume $X$ to be a space of functions on a one-dimensional interval, which we discretize with a mesh of $c$ equidistant points. The results provided in \crefsection:simulation_results include two-dimensional examples as well. Discretizing $F$ accordingly, here using a finite difference method, we obtain from \eqrefeq:first-order_system_discrete_in_time the fully discrete evolution

Each element in the vector ${\mathbf{z}}_{k}$ represents an approximation of the true state $z$ at one of the spatial grid points and time $t=k\,\Delta t$.

A neural-network predictor takes the current state vector ${\mathbf{z}}_{k}$ (or a reconstruction thereof) and predicts the solution at the next time step ${\mathbf{z}}_{k+1}$. The estimator network is trained over simulation data or acquired real-life data. In most cases, the state ${\mathbf{z}}_{k}$ cannot be fully measured by an array of sensors. Therefore, the data provided by the sensors available represents an incomplete state vector and the full state vector is reconstructed first before being fed into the predictor network.

The reconstruction of an incomplete state reading works as follows. A sensor arrangement is encoded in a $c$-dimensional vector $\omega$. Each element of $\omega$ indicates the weight of the presence of a sensor at the respective vertex. Weights less than 0.5 will be treated as no sensor present. The function $f_{r}\colon\R^{c}\times\R^{c}\to\R^{c}$ provides an approximate reconstruction based on an interpolation of neighboring sensors measurements, \Creffig-recons. That is, we have ${\mathbf{z}}\approx f_{r}({\mathbf{z}},\omega)$. A cubic spline interpolation is used in the one-dimensional simulations. This interpolation is a suitable choice for PDEs with second-order spatial derivatives.

The function $f_{n}\colon\R^{c}\to\R^{c}$ provides a neural-network prediction ${\mathbf{Z}}_{k+1}$ of the state ${\mathbf{z}}_{k+1}$ using reconstructed data; that is

Different types of neural networks lead to different formulations for the predictor function $f_{n}$. Examples are given in \crefsection:simulation_results.

In this section, we consider an optimization of the sensor configuration encoded in the weight vector $\omega$. Clearly, from the prediction point of view, it is beneficial to have sensors everywhere. In this case, we would have ${\mathbf{z}}_{k}=f_{r}({\mathbf{z}}_{k},\omega)$ since no interpolation is necessary. On the other hand, we associate a cost with each sensor placed. Our goal is to obtain a compromise between the number of sensors and the accuracy of the predictions ${\mathbf{Z}}_{k+1}$ of the state ${\mathbf{z}}_{k+1}$, based on the partially observed previous state ${\mathbf{z}}_{k}$. We express this goal in terms of the cost function {multline} J(ω) = α  ∑_i=1^c[ω]_i
+ ∑_k=0^K-1 ∑_i=1^c \paren[big][]f_n(f_r(z_k,ω)) - z_k+1_i^2   h_i   Δt . Here $K$ denotes the number of grid points in time. Notice that the second term in \eqrefeq:cost_function represents a discretized version of the $L^{2}$ mean squared error. The sequence $h_{i}$ depends on the spatial mesh size and the quadrate scheme used; for instance, $h_{i}=h/2,h,\cdots,h,h/2$ implements the trapezoidal rule with mesh size $h$. A similar formula holds in two spatial dimensions.

The derivative of $J(\omega)$ with respect to (\wrt) the sensor weight vector $\omega$ is {multline} [J’(ω)]_i = α
+ 2∑_k=0^K-1 ∑_i=1^c \paren[big][]\paren[big]()f_n(f_r(z_k,ω)) - z_k+1^\transpZ_k+1’_i   h_i   Δt . Here ${\mathbf{Z}}_{k+1}^{\prime}$ is the derivative of ${\mathbf{Z}}_{k+1}=f_{n}(f_{r}({\mathbf{z}}_{k},\omega))$ \wrt $\omega$. From the chain rule, we obtain

where the right hand side is a product of two $c\times c$ matrices. The first factor represents the derivative of the estimator network’s output \wrt its input vector, while the second factor is the derivative of the outcome of the spline interpolation \wrt $\omega$. The latter derivative is zero except when there are components of $\omega$ right at the threshold $\omega=0.5$. In this case, the derivative becomes a Dirac function. Simulation results, however, indicate that the first term in \eqrefeq:chain_rule is sufficient to provide search directions suitable to reduce the cost function. Therefore, we use the approximation

to evaluate the derivative $J^{\prime}(\omega)$ in \eqrefeq:derivative_of_cost_function. Its transpose, the approximate gradient $\nabla J(\omega)$, is then used in a trust-region optimization algorithm to find an improved sensor configuration vector $\omega$.

In this section, we provide numerical results for a range of different PDE models. Our goal is to demonstrate that, in each case, optimized sensor configurations and properly trained predictor networks $f_{n}$ can obtain sufficiently accurate predictions of the state vector. To this end, we show that the cost $J(\omega)$ for an optimized weight $\omega$ is significantly lower than the cost for the all-sensor case $[\omega]_{i}=1$, $i=1,2,\ldots,c$, for which the lower bound $J(\omega)\geq\alpha\,c$ is valid. The overall procedure is described in \crefalgorithm:deep_learning_and_shape_optimization.

For each example, we build an estimator neural network $f_{n}$ using the deep learning package Keras (version 2.4.0). We employ a linear activation function and the mean squared error (MSE) loss function in the networks and use the Adam optimizer for training. The Adam optimizer and the MSE loss function yielded better performance compared to the rest of available optimization algorithms and loss functions. For PDE models in one space dimension, the network layout consists of two dense layers, where each layer has as many neurons as grid nodes. For PDE models in two space dimensions, a convolutional neural network with one layer, one filter and kernel size $3\times 3$ and zero padding is considered. In each case, the gradient of the network is calculated using the method backend in Keras.

In the implementation of the reconstructor function $f_{r}$, the sensor data is reconstructed using the Python package interpolate. For PDE models in one space dimension, CubicSpline is used, whereas interp2d is chosen in two space dimensions.

Simulations are presented for a variety of PDE models in the following subsections. In figures and formulas, $u_{r}$ denotes the true (simulated) solution, which corresponds to ${\mathbf{z}}$ in the notation of \crefsection:estimator_design,section:shape_optimization. Moreover, $u_{p}$ denotes the predicted solution corresponding to ${\mathbf{Z}}$ according to \eqrefeq:reconstruction_followed_by_prediction.

Optimal prediction for distributed parameter systems was discussed in this paper. The optimal predictions involves two steps. In the first step, a neural network predictor is designed. In the second step, sensor shapes are optimized using the gradient of the network. Computer simulations are conducted for several models including heat and wave equation in one and two space dimensions. The results show a significant reduction in cost of prediction as well as improvement in the predictor performance. The computer codes are accessible in the repository \urlhttps://github.com/TUCMath/optimal-predictor.