Toward real-time optimization through model reduction and model discrepancy sensitivities

TL;DR

This paper introduces a method combining hyper-differential sensitivity analysis with model discrepancy to enhance the reliability of reduced-order model constrained optimization, enabling near real-time solutions for complex systems.

Contribution

It proposes a novel HDSA-MD approach that improves ROM-based optimization accuracy with minimal full-order model evaluations during online phases.

Findings

HDSA-MD significantly improves optimization solutions.

Only one full-order model evaluation needed online.

Effective in atmospheric contaminant control and wildfire estimation.

Abstract

Optimization problems arise in a range of scenarios, from optimal control to model parameter estimation. In many applications, such as the development of digital twins, it is essential to solve these optimization problems within wall-clock-time limitations. However, this is often unattainable for complex systems, such as those modeled by nonlinear partial differential equations. One strategy for mitigating this issue is to construct a reduced-order model (ROM) that enables more rapid optimization. In particular, the use of nonintrusive ROMs -- those that do not require access to the full-order model at evaluation time -- is popular because they facilitate optimization solutions can be computed within the wall-clock-time requirements. However, the optimization solution will be unreliable if the iterates move outside the ROM training data. This article proposes the use of hyper-differential sensitivity analysis with respect to model discrepancy (HDSA-MD) as a computationally efficient tool to augment ROM-constrained optimization and improve its reliability. The proposed approach consists of two phases: (i) an offline phase where several full-order model evaluations are computed to train the ROM, and (ii) an online phase where a ROM-constrained optimization problem is solved, $N=\mathcal{O}(1)$ full-order model evaluations are computed, and HDSA-MD is used to enhance the optimization solution using the full-order model data. Numerical results are demonstrated for two examples, atmospheric contaminant control and wildfire ignition location estimation, in which a ROM is trained offline using inaccurate atmospheric data. The HDSA-MD update yields a significant improvement in the ROM-constrained optimization solution using only one full-order model evaluation online with corrected atmospheric data.

Tables

Table 1. Table 1: FOM misfit (normalized by ‖ 𝒀 ‖ F 2 \|\boldsymbol{Y}\|^{2}_{F} ) from running WRF-SFIRE at the ROMCO solution 𝒛 ~ \tilde{{\boldsymbol{z}}} and the posterior mean ignition location 𝒛 ¯ = 𝒛 ~ − 𝑯 − 1 ​ 𝑩 ​ 𝜽 ¯ \overline{{\boldsymbol{z}}}=\tilde{{\boldsymbol{z}}}-{\boldsymbol{H}}^{-1}{\boldsymbol{B}}\overline{{\boldsymbol{\theta}}} .

Test Case	Relative misfit at $\widetilde{𝒛}$	Relative misfit at $\overline{𝒛}$	Reduction
1	$1.25\times {10}^{-3}$	$7.98\times {10}^{-4}$	36%
2	$1.90\times {10}^{-1}$	$8.23\times {10}^{-2}$	57%
3	$1.74\times {10}^{-3}$	$4.51\times {10}^{-4}$	74%
4	$3.67\times {10}^{-2}$	$1.22\times {10}^{-2}$	67%
5	$4.13\times {10}^{-2}$	$1.71\times {10}^{-2}$	59%
Average	$5.42\times {10}^{-2}$	$2.26\times {10}^{-2}$	58%

Table 2. Table 2: Error (in meters) of the ROMCO solution 𝒛 ~ \tilde{{\boldsymbol{z}}} and the posterior mean ignition location 𝒛 ¯ \overline{{\boldsymbol{z}}} ; 𝒛 † {\boldsymbol{z}}^{\dagger} denotes the “ground truth” ignition location used in the test set simulation.

Test Case	$\Vert \widetilde{𝒛}-{𝒛}^{†}\Vert$	$\Vert \overline{𝒛}-{𝒛}^{†}\Vert$	Reduction
1	$2.81\times {10}^2$	$2.31\times {10}^2$	18%
2	$2.14\times {10}^3$	$1.41\times {10}^3$	34%
3	$3.26\times {10}^2$	$1.33\times {10}^2$	59%
4	$9.99\times {10}^2$	$5.45\times {10}^2$	45%
5	$1.08\times {10}^3$	$6.60\times {10}^2$	39%
Average	$9.67\times {10}^2$	$5.96\times {10}^2$	38%

Equations

$$\underset{\mathbf{z}\in {\mathbb{R}}^{n_z}}{\min }J(\mathbf{S}(t;\mathbf{z},\mathbf{\xi }),\mathbf{z})$$

$$\text{where }\mathbf{S}(t;\mathbf{z},\mathbf{\xi })\text{ solves}$$

$$\begin{array}{rl}& \dot{\mathbf{u}}(t)=\mathbf{f}(\mathbf{u}(t),\mathbf{z},\mathbf{\xi })t\in [0,T]\\ & \mathbf{u}(0)=\mathring{\mathbf{u}}.\end{array}$$

$$\mathbf{u}(t)\approx \mathbf{V}\hat{\mathbf{u}}(t)=\sum _{i=1}^{n_{\hat{u}}}{\hat{u}}_i(t){\mathbf{v}}_i.$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^{n_z}}{\min }J(\mathbf{V}\hat{\mathbf{S}}(t;\mathbf{z}),\mathbf{z})$$

$$\text{where }\hat{\mathbf{S}}(t;\mathbf{z})\text{ solves}$$

$$\begin{array}{rl}& \dot{\hat{\mathbf{u}}}(t)=\hat{\mathbf{f}}(\hat{\mathbf{u}}(t),\mathbf{z})t\in [0,T]\\ & \hat{\mathbf{u}}(0)={\mathbf{V}}^{\mathsf{T}}\mathring{\mathbf{u}}.\end{array}$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^{n_z}}{\min }J(\mathbf{V}\hat{\mathbf{S}}(t;\mathbf{z})+\mathbf{\delta }(t;\mathbf{z},\mathbf{\theta }),\mathbf{z})$$

$$\text{where }\hat{\mathbf{S}}(t;\mathbf{z})\text{ solves}$$

$$\begin{array}{rl}& \dot{\hat{\mathbf{u}}}(t)=\hat{\mathbf{f}}(\hat{\mathbf{u}}(t),\mathbf{z})t\in [0,T]\\ & \hat{\mathbf{u}}(0)={\mathbf{V}}^{\mathsf{T}}\mathring{\mathbf{u}}.\end{array}$$

$$\displaystyle{\boldsymbol{\delta}}(t;{\boldsymbol{z}}_{\ell},{\boldsymbol{\theta}})\approx{\boldsymbol{S}}(t;{\boldsymbol{z}}_{\ell},{\boldsymbol{\xi}}^{\dagger})-{\boldsymbol{V}}\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}}_{\ell}),\qquad\ell=1,2,\dots,N.\$$

$$\mathbf{\delta }(t;\mathbf{z},\mathbf{\theta })={\mathbf{u}}_0(t;\mathbf{\theta })+\mathbf{L}(t;\mathbf{\theta }){\mathbf{M}}_{\mathbf{z}}\mathbf{z}.$$

$${\nabla }_{\mathbf{z}}\mathcal{J}(\tilde{\mathbf{z}},{\mathbf{\theta }}_0)=0,$$

$$\mathcal{J}(\mathbf{z},\mathbf{\theta })=J(\mathbf{V}\hat{\mathbf{S}}(t;\mathbf{z})+\mathbf{\delta }(t;\mathbf{z},\mathbf{\theta }),\mathbf{z})$$

$${\nabla }_{\mathbf{z}}\mathcal{J}({\mathbf{z}}^{\star }(\mathbf{\theta }),\mathbf{\theta })=0$$

$${\nabla }_{\mathbf{z}}{\mathbf{z}}^{\star }({\mathbf{\theta }}_0)=-{\mathbf{H}}^{-1}\mathbf{B}$$

$${\mathbf{z}}^{\star }(\mathbf{\theta })\approx {\mathbf{z}}^{\star }({\mathbf{\theta }}_0)+{\nabla }_{\mathbf{\theta }}{\mathbf{z}}^{\star }({\mathbf{\theta }}_0)\mathbf{\theta }=\tilde{\mathbf{z}}-{\mathbf{H}}^{-1}\mathbf{B}\mathbf{\theta }$$

$$\frac{\partial u_1}{\partial t}={\kappa }_1\Delta u_1-\mathbf{v}\cdot \nabla u_1-\rho u_1{u}_2$$

$$\frac{\partial u_2}{\partial t}={\kappa }_2\Delta u_2-\mathbf{v}\cdot \nabla u_2-\rho u_1{u}_2+s$$

$$\mathbf{n}\cdot \nabla u_1=\mathbf{n}\cdot \nabla u_2=0$$

$$u_1(\mathbf{x},0)={\mathring{u}}_1(\mathbf{x}),u_2(\mathbf{x},0)=0$$

$$\mathbf{M}{\dot{\mathbf{u}}}_1(t)$$

$$\mathbf{M}{\dot{\mathbf{u}}}_2(t)$$

$${\mathbf{u}}_1(0)$$

$$s(\mathbf{x},t)=\sum _{i=1}^{n_q}{q}_i(t;\mathbf{z}{)}^2{\varphi }_i(\mathbf{x}),{\varphi }_i(\mathbf{x})=50\exp \left(-1000\Vert \mathbf{x}-{\bar{\mathbf{x}}}_i{\Vert }^2\right),$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^{n_z}}{\min }\frac{1}{2}{\int }_0^T\Vert {\mathbf{u}}_1(t)\odot \mathbf{\psi }{\Vert }_{\mathbf{M}}^2+\gamma \Vert \mathbf{q}(t;\mathbf{z})\odot \mathbf{q}(t;\mathbf{z}){\Vert }^{2}dt$$

$$\text{where}(\mathbf{u}(t),\mathbf{q}(t;\mathbf{z}))\text{ jointly solves lx@crefcreftype refnumeq:adr-ode},$$

$$\displaystyle{\boldsymbol{v}}({\boldsymbol{x}};{\boldsymbol{\xi}})=\left[\begin{array}[]{c}4\xi_{1}\cos(\xi_{4}x_{1}/L)-2\xi_{2}\sin(\xi_{3}\pi x_{1})\sin(\xi_{4}x_{1}/L)\\ 4\xi_{1}\sin(\xi_{4}x_{1}/L)-2\xi_{2}\sin(\xi_{3}\pi x_{1})\cos(\xi_{4}x_{1}/L)\end{array}\right].$$

$${\mathbf{u}}_1(t)\approx {\mathbf{V}}_1{\hat{\mathbf{u}}}_1(t),{\mathbf{u}}_2(t)\approx {\mathbf{V}}_2{\hat{\mathbf{u}}}_2(t),$$

$${\dot{\hat{\mathbf{u}}}}_1(t)$$

$${\dot{\hat{\mathbf{u}}}}_2(t)$$

$$\begin{array}{rl}{\hat{\mathbf{A}}}_k& ={\mathbf{V}}_k^{\mathsf{T}}{\mathbf{A}}_k{\mathbf{V}}_k\in {\mathbb{R}}^{r_k\times r_k},\\ {\hat{\mathbf{R}}}_k& ={\mathbf{V}}_k^{\mathsf{T}}{\mathbf{R}}_k({\mathbf{V}}_1\otimes {\mathbf{V}}_2)\in {\mathbb{R}}^{r_k\times r_1{r}_2},\end{array}k=1,2,\text{and}\hat{\mathbf{\Phi }}={\mathbf{V}}^{\mathsf{T}}\mathbf{M}\mathbf{\Phi }\in {\mathbb{R}}^{r_2\times n_q}.$$

$$\underset{\mathbf{z}\in {\mathbb{R}}^{n_z}}{\min }\frac{1}{2}{\int }_0^T\Vert {\mathbf{V}}_1{\hat{\mathbf{u}}}_1(t)\odot \mathbf{\psi }{\Vert }_{\mathbf{M}}^2+\gamma \Vert \mathbf{q}(t;\mathbf{z})\odot \mathbf{q}(t;\mathbf{z}){\Vert }^{2}dt$$

$$\text{where}(\hat{\mathbf{u}}(t),\mathbf{q}(t;\mathbf{z}))\text{ jointly solves lx@crefcreftype refnumeq:adr-rom},$$

$$\mathcal{J}(\mathbf{z},{\mathbf{\theta }}_0)\approx \mathcal{J}(\tilde{\mathbf{z}},{\mathbf{\theta }}_0)+\frac{1}{2}(\mathbf{z}-\tilde{\mathbf{z}}{)}^T\mathbf{H}(\mathbf{z}-\tilde{\mathbf{z}})$$

$$\tilde{\mathbf{z}}-\mathbf{P}{\mathbf{H}}^{-1}\mathbf{B}\mathbf{\theta }.$$

$$\frac{\partial u}{\partial t}+\Psi (\mathbf{v})\Vert \nabla u\Vert =0,$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Advancements in modeling and computing capabilities have opened up the possibility of creating digital twins for complex systems. This can impact a range of applications including additive manufacturing, fusion energy, aerospace systems, and wildfire containment. Stated briefly, a digital twin is an goal-oriented coupling of physical and virtual assets with a bidirectional information flow between the assets [28]. Central to the digital twin concept is the ability to adapt and inform decision-making as data is collected in real-time—which is defined relative to the timescale of the data collection and decisions. This requires solving optimization problems constrained by numerical models to (i) calibrate the digital twin to physical system data, and (ii) use the digital twin to inform decisions, such as actions to control the physical system. However, the computational cost of model evaluations makes real-time optimization unattainable in many applications. This article seeks to enable real-time optimization for systems modeled by partial differential equations (PDEs) through a novel combination of non-intrusive model reduction and hyper-differential sensitivity analysis with respect to model discrepancy.

Model reduction seeks to alleviate the computational burden of solving a high-dimensional system, referred to as the full-order model (FOM), by identifying and exploiting structure to build a low-dimensional system, referred to as the reduced-order model (ROM). Classical model reduction is intrusive in the sense that it learns low-dimensional structure from training data but also requires access to the FOM to construct and/or solve the ROM. This is unsuitable for many applications requiring real-time computation because the cost of FOM access inside the optimization loop remains prohibitively expensive. Non-intrusive model reduction methods seek to overcome this limitation by learning a low-dimensional system exclusively from training data so that additional access to the FOM is not needed. This ensures that the computational cost of reduced-order model evaluations is independent of the FOM dimension. However, this gain in computational efficiency generally results in a loss in accuracy, especially when predicting outside of the regime represented by the training data. This presents a challenge for optimization problems constrained by the ROM, as the optimizer’s iterations may drift away from the training data regime, leading to an unreliable solution.

To improve the reliability of ROM-constrained optimization, we propose to augment it with hyper-differential sensitivity analysis with respect to model discrepancy, herein referred to as HDSA-MD. The approach originated in [15, 16] to facilitate optimization for problems where a high-fidelity model can only be evaluated a small number of times. Specifically, these articles consider optimization constrained by a low-fidelity model, derived from physics simplifications such as omission of small-magnitude forces or linearization, and use $N=\mathcal{O}(1)$ evaluations of a high-fidelity PDE model to update the low-fidelity optimization solution. The benefits of HDSA-MD are that high-fidelity derivative information is not required and the evaluations can be performed outside of the optimization loop. This article combines HDSA-MD and non-intrusive reduced-order modeling to enable real-time optimization. Computational efficiency is achieved using ROM-constrained optimization that is independent of the FOM, and reliability is improved via the HDSA-MD optimal solution update.

We demonstrate our approach by training an operator inference ROM to enable real-time atmospheric contaminant control and by training a neural operator ROM to enable real-time wildfire ignition source estimation. In both examples, we require that the ROMs are trained “offline” with inaccurate estimates of the atmospheric conditions. These pre-trained ROMs are insufficient since the state of the atmosphere at decision time is generally different than in the training data. However, by augmenting the ROM-constrained optimization with HDSA-MD, we show that using as little as one FOM evaluation with accurate atmospheric data is sufficient to enhance the ROM-constrained optimization.

There is a vast literature on the use of ROMs to reduce the computational cost of PDE-constrained optimization. Goal-oriented model reduction in the service of optimization is introduced in [7]. To improve ROM accuracy, [43, 34, 44, 13] propose various approaches that construct a sequence of optimization problems and iteratively enrich the ROM. Specialized ROMs for Hessian approximation are considered in [18]. ROMs are used to enable shape optimization in [2, 3] and risk-averse optimization in [19, 45]. Recent work has also explored the use of derivative information to train neural operator ROMs in the service of solving optimization problems [30, 24]. This article differs from previous works in that we focus on real-time applications where a non-intrusive ROM is trained offline in preparation for ROM-constrained optimization online, yet a small number (as little as 1) of FOM evaluations are feasible online to enhance the optimization solution using HDSA-MD.

The article is organized as follows. Section 2 presents the optimization problem under consideration and the corresponding ROM-constrained optimization problem derived from it. We present our HDSA-MD formulation for ROM-constrained optimization in Section 3. Section 4 and Section 5 demonstrate our approach for ROM-constrained optimization arising in control of atmospheric contaminants and estimation of wildfire ignition, respectively. Concluding remarks are made in Section 6.

2 Optimization formulation

Consider FOM-constrained optimization problems of the form

[TABLE]

Time is denoted by $t\in[0,T]$, $T>0$ is the final time, ${\boldsymbol{u}}(t)\in{\mathbb{R}}^{n_{u}}$ is the state, ${\boldsymbol{z}}\in{\mathbb{R}}^{n_{z}}$ corresponds to control, design, or inversion parameters, ${\boldsymbol{\xi}}\in{\mathbb{R}}^{n_{\xi}}$ is a vector of model parameters, ${\boldsymbol{f}}:{\mathbb{R}}^{n_{u}}\times{\mathbb{R}}^{n_{z}}\times{\mathbb{R}}^{n_{\xi}}\to{\mathbb{R}}^{n_{u}}$ defines the system dynamics, $\mathring{{\boldsymbol{u}}}\in{\mathbb{R}}^{n_{u}}$ is the initial state, and ${\boldsymbol{S}}(t;{\boldsymbol{z}},{\boldsymbol{\xi}})$ is the solution operator that maps $({\boldsymbol{z}},{\boldsymbol{\xi}})$ to the solution of (2.1b). We assume that ${\boldsymbol{S}}(t;{\boldsymbol{z}},{\boldsymbol{\xi}})$ is uniquely defined for each $({\boldsymbol{z}},{\boldsymbol{\xi}})$. The objective function $J:L^{2}([0,T])\times{\mathbb{R}}^{n_{z}}\to{\mathbb{R}}$ depends on the solution operator ${\boldsymbol{S}}$, and hence each objective function evaluation requires solving (2.1b). We seek to solve FOMCO with fixed model parameters ${\boldsymbol{\xi}}$; however, the model parameters are uncertain, so variability of ${\boldsymbol{\xi}}$ plays an important role in our analysis. In our numerical results, ${\boldsymbol{\xi}}$ corresponds to atmospheric winds, which are only known within a short time window before the optimization solution is required.

The FOMCO problem arises in many applications where a PDE is discretized in space, yielding the state vector ${\boldsymbol{u}}(t)\in{\mathbb{R}}^{n_{u}}$. In this case, $n_{u}$ is the number of nodes in the spatial discretization, making eq. 2.1b a high-dimensional system. Throughout the article we will refer to ${\boldsymbol{z}}\in{\mathbb{R}}^{n_{z}}$ as the decision variable, noting that it may correspond to control, design, or inversion parameters. We formulate FOMCO with the state written as a function of time, i.e., ${\boldsymbol{u}}(t)$, but the decision variables ${\boldsymbol{z}}$ and model parameters ${\boldsymbol{\xi}}$ as vectors in Euclidean space. In general, we may consider ${\boldsymbol{z}}$ and ${\boldsymbol{\xi}}$ that are functions of space and/or time, but we assume that they are discretized in FOMCO to simplify the presentation. Where it is relevant, we will use weighted inner products to bring information from the function spaces and discretization into our analysis.

Assuming that FOMCO cannot be solved within wall-clock-time requirements, we consider a process in which FOM evaluations are performed offline in preparation for online computation to approximate the solution of FOMCO. In the offline phase, we generate data by evaluating the FOM for $M$ decision variables instances ${\boldsymbol{z}}_{k}\in{\mathbb{R}}^{n_{z}}$, $k=1,2,\dots,M$, to produce a set of state trajectories ${\boldsymbol{u}}_{k}(t)$, $k=1,2,\dots,M$, satisfying (2.1b). From these trajectories, a low-dimensional linear subspace, $\text{span}\{{\boldsymbol{v}}_{1},{\boldsymbol{v}}_{2},\dots,{\boldsymbol{v}}_{n_{\hat{u}}}\}$, is determined, where ${\boldsymbol{v}}_{i}\in{\mathbb{R}}^{n_{u}}$, $i=1,2,\dots,n_{\hat{u}}$, are basis vectors computed from FOM data. We let ${\boldsymbol{V}}=\begin{pmatrix}{\boldsymbol{v}}_{1}&{\boldsymbol{v}}_{2}&\cdots&{\boldsymbol{v}}_{n_{\hat{u}}}\end{pmatrix}\in{\mathbb{R}}^{n_{u}\times n_{\hat{u}}}$ denote the basis matrix and represent the state ${\boldsymbol{u}}(t)$ by the time-dependent coordinates $\widehat{{\boldsymbol{u}}}(t)=(\hat{u}_{1}(t),\ldots,\hat{u}_{n_{\hat{u}}}(t))\in\mathbb{R}^{n_{\hat{u}}}$ as

[TABLE]

This reduces the state degrees-of-freedom from $n_{u}$ to $n_{\hat{u}}$, where $n_{\hat{u}}\ll n_{u}$. The full dynamics ${\boldsymbol{f}}:{\mathbb{R}}^{n_{u}}\times{\mathbb{R}}^{n_{z}}\times{\mathbb{R}}^{n_{\xi}}\to{\mathbb{R}}^{n_{u}}$ are approximated by reduced dynamics $\widehat{{\boldsymbol{f}}}:{\mathbb{R}}^{n_{\hat{u}}}\times{\mathbb{R}}^{n_{z}}\times{\mathbb{R}}^{n_{\xi}}\to{\mathbb{R}}^{n_{\hat{u}}}$. There are a plurality of ways to construct $\widehat{{\boldsymbol{f}}}$. For instance, Galerkin projection defines $\widehat{{\boldsymbol{f}}}(\widehat{{\boldsymbol{u}}},{\boldsymbol{z}},{\boldsymbol{\xi}})={\boldsymbol{V}}^{\mathsf{T}}{\boldsymbol{f}}({\boldsymbol{V}}\widehat{{\boldsymbol{u}}},{\boldsymbol{z}},{\boldsymbol{\xi}})$. However, this is intrusive in the sense that evaluating $\widehat{{\boldsymbol{f}}}$ requires access to ${\boldsymbol{f}}$. This article focuses on non-intrusive methods—those which only require data output from FOM solves—to construct $\widehat{{\boldsymbol{f}}}$. In particular, we demonstrate our framework using operator inference [32] and neural operator surrogates [20], which are detailed in Section 4 and Section 5, respectively.

Due to the complexity of learning a non-intrusive ROM from limited FOM data, we consider the case where the model parameters ${\boldsymbol{\xi}}$ are fixed to a nominal value $\tilde{{\boldsymbol{\xi}}}$ for FOM data generation and ROM training. We suppress the dependence of $\widehat{{\boldsymbol{f}}}$ on ${\boldsymbol{\xi}}$ with the implication that $\widehat{{\boldsymbol{f}}}$ is trained using data with ${\boldsymbol{\xi}}=\tilde{{\boldsymbol{\xi}}}$. We will account for uncertainty in ${\boldsymbol{\xi}}$ through our proposed HDSA-MD approach.

Given a basis matrix ${\boldsymbol{V}}$ and approximate dynamics $\widehat{{\boldsymbol{f}}}$, we consider the ROM-constrained optimization problem

[TABLE]

Thanks to the dimension reduction $n_{\hat{u}}\ll n_{u}$, the ROM solution operator $\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}})$ can be evaluated much faster than the FOM solution operator ${\boldsymbol{S}}(t;{\boldsymbol{z}},{\boldsymbol{\xi}})$ and thereby enable real-time performance.

Let ${\tilde{{\boldsymbol{z}}}}$ denote the ROMCO minimizer that is computed online, but that is insufficient due to the discrepancy between the FOM and ROM. This discrepancy stems from three sources: (i) restriction to the basis matrix ${\boldsymbol{V}}$; (ii) failing to learn the complete FOM dynamics; and (iii) misspecification of ${\boldsymbol{\xi}}$ in training. Our goal in this article is to improve the ROMCO solution ${\tilde{{\boldsymbol{z}}}}$ in real-time using $N$ FOM solves. In particular, we focus on the case where only $N=\mathcal{O}(1)$ FOM evaluations are achievable within wall-clock-time limitations. We assume that the “true” model parameters ${\boldsymbol{\xi}}={\boldsymbol{\xi}}^{\dagger}$ are available online so that the $N$ FOM solves use the model parameters ${\boldsymbol{\xi}}^{\dagger}$ rather than the model parameters $\tilde{{\boldsymbol{\xi}}}$ used to generate data for ROM training. The goal is to update ${\tilde{{\boldsymbol{z}}}}$ to better approximate the FOMCO solution when ${\boldsymbol{\xi}}={\boldsymbol{\xi}}^{\dagger}$.

3 HDSA with respect to model discrepancy

Consider the model discrepancy parameterized optimization problem

[TABLE]

ROMCO-MD corresponds to endowing ROMCO with a model discrepancy operator ${\boldsymbol{\delta}}(t;{\boldsymbol{z}},{\boldsymbol{\theta}})$, parameterized by ${\boldsymbol{\theta}}\in{\mathbb{R}}^{n_{\theta}}$, which seeks to account for the discrepancy between FOM and ROM, i.e., ${\boldsymbol{S}}(t;{\boldsymbol{z}},{\boldsymbol{\xi}}^{\dagger})-{\boldsymbol{V}}\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}})$. In what follows, we demonstrate how to use this model discrepancy endowment to enhance the ROMCO solution.

3.1 Model discrepancy calibration

Assume that the FOM is evaluated $N$ times to provide discrepancy data $\{{\boldsymbol{z}}_{\ell},{\boldsymbol{S}}(t;{\boldsymbol{z}}_{\ell},{\boldsymbol{\xi}}^{\dagger})-{\boldsymbol{V}}\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}}_{\ell})\}_{\ell=1}^{N}$. We seek to parameterize ${\boldsymbol{\delta}}$ and calibrate its parameters ${\boldsymbol{\theta}}$ so that

[TABLE]

In general, the discrepancy between FOM and ROM depends nonlinearly on ${\boldsymbol{z}}$. However, we cannot expect to learn this nonlinearity since the training dataset only contains $N=\mathcal{O}(1)$ FOM evaluations. Rather, the model discrepancy operator is parameterized as affine in ${\boldsymbol{z}}$:

[TABLE]

Here, ${\boldsymbol{u}}_{0}:[0,T]\times{\mathbb{R}}^{n_{\theta}}\to{\mathbb{R}}^{n_{u}}$ and $\boldsymbol{L}:[0,T]\times{\mathbb{R}}^{n_{\theta}}\to{\mathbb{R}}^{n_{u}\times n_{z}}$ are the intercept and linear terms in the affine operator. The mass matrix ${\boldsymbol{M}}_{\boldsymbol{z}}\in{\mathbb{R}}^{n_{z}\times n_{z}}$, whose $(i,j)$ entry is the inner product of the $i^{th}$ and $j^{th}$ decision-variable basis function111The mass matrix ${\boldsymbol{M}}_{\boldsymbol{z}}$ is the identity matrix if ${\boldsymbol{z}}$ does not arise from discretization of a function., ensures that the discrepancy parameterization respects the structure of the decision variables if they arise from spatial and/or temporal discretization. The intercept ${\boldsymbol{u}}_{0}$ and linear operator $\boldsymbol{L}$ are parameterized as linear functions of ${\boldsymbol{\theta}}$ and hence ${\boldsymbol{\delta}}(t;{\boldsymbol{z}},{\boldsymbol{\theta}})$ depends linearly on ${\boldsymbol{\theta}}$.

To calibrate the model discrepancy parameters, we first observe that the dimension of ${\boldsymbol{\theta}}\in{\mathbb{R}}^{n_{\theta}}$ scales with the product of the state and decision variable dimensions. If the number of FOM evaluations is less than the decision variable dimension, i.e., $N<n_{z}$, then there exists infinitely many ${\boldsymbol{\theta}}$’s that fit the data (3.2). Accordingly, we adopt a Bayesian calibration paradigm so that prior knowledge can constrain the calibration and uncertainty can be accounted for through the posterior. Bayesian calibration is challenging because ${\boldsymbol{\delta}}$ is an operator mapping from the decision variable space to the state space. We summarize the approach below and refer the reader to [16, 17] for additional details.

Prior discrepancy

The prior for ${\boldsymbol{\theta}}$ should be defined on the space of affine operators and discretized to ensure that it respects the pre-discretization structures. Following [16, 17], we consider a mean-zero Gaussian prior for the discrepancy parameters ${\boldsymbol{\theta}}$. Recall that, for a general a mean-zero Gaussian random vector $\boldsymbol{X}$, its probability density function is proportional to $\exp\left(-\frac{1}{2}\|\boldsymbol{x}\|_{{\boldsymbol{W}}_{\boldsymbol{x}}}\right)$, where ${\boldsymbol{W}}_{\boldsymbol{x}}$ is the precision matrix (the inverse of the covariance matrix). We define the prior covariance for the discrepancy parameters ${\boldsymbol{\theta}}$ by introducing a weighted norm $\|\cdot\|_{{\boldsymbol{W}}_{\boldsymbol{\theta}}}$ and take the weighting matrix ${\boldsymbol{W}}_{\boldsymbol{\theta}}$ as the prior precision.

To derive ${\boldsymbol{W}}_{\boldsymbol{\theta}}$, we define mean-zero Gaussian priors on the state and decision variables spaces, resulting in precision matrices ${\boldsymbol{W}}_{\boldsymbol{u}}$ and ${\boldsymbol{W}}_{\boldsymbol{z}}$ that define inner products for ${\boldsymbol{u}}$ and ${\boldsymbol{z}}$, respectively. Considering an inner product on the space of affine operators, we can derive an expression for ${\boldsymbol{W}}_{\boldsymbol{\theta}}$ that depends on ${\boldsymbol{W}}_{\boldsymbol{u}}$ and ${\boldsymbol{W}}_{\boldsymbol{z}}$. These precision matrices are defined using theory and experience from the analysis of infinite-dimensional Bayesian inverse problems [41, 6]. See [17] for detailed analysis of the prior specification.

Posterior discrepancy

The data likelihood function is defined by assuming a mean-zero Gaussian noise model. Since the data consists of full spatio-temporal observations of the state discrepancy, we define the noise covariance as a scalar multiple of the identity operator in the state space. This corresponds to a Gaussian noise model with covariance $\alpha_{d}{\boldsymbol{M}}_{\boldsymbol{u}}^{-1}$, where $\alpha_{d}>0$ scales the variance and ${\boldsymbol{M}}_{\boldsymbol{u}}$ defines the state-space inner product222The matrix ${\boldsymbol{M}}_{\boldsymbol{u}}$ defines an inner product that integrates over both time and space. Hence, ${\boldsymbol{M}}_{\boldsymbol{u}}$ does not correspond to the mass matrix in the spatial discretization, but rather is a Kronecker product of a temporal weighting matrix and the spatial discretization mass matrix.. In Bayesian inverse problems, the noise level is determined based on characteristics of the data-collection process. In our context, the data is not noisy as is typical with sensor measurements. Instead, data arises from simulations, and therefore the level of noise corresponds to linearization error from fitting an affine model to data from a nonlinear operator. We interpret $\alpha_{d}$ as a linearization error magnitude and hence specify it relative to the magnitude of the discrepancy data being fit.

Since we have assumed Gaussian prior and noise models, and parameterized ${\boldsymbol{\delta}}$ as a linear function of ${\boldsymbol{\theta}}$, the posterior distribution for ${\boldsymbol{\theta}}$ is Gaussian, with its mean and covariance matrix available in closed form. To simplify the exposition, we omit explicit expressions for the posterior mean and samples, but refer the reader to [16], which demonstrates efficient posterior computation.

3.2 Optimal solution updating

Our goal is to approximate the FOMCO solution online using the ROMCO minimizer ${\tilde{{\boldsymbol{z}}}}$ and $N=\mathcal{O}(1)$ FOM evaluations. To analyze how the model discrepancy influences the optimization problem, observe that if ${\boldsymbol{\theta}}_{0}=\boldsymbol{0}$, then ${\boldsymbol{\delta}}(t;{\boldsymbol{z}},{\boldsymbol{\theta}}_{0})\equiv 0$ and the ROMCO minimizer ${\tilde{{\boldsymbol{z}}}}$ satisfies the first-order optimality condition

[TABLE]

where

[TABLE]

is the objective function of the ROMCO-MD problem (which coincides with the ROMCO objective function when ${\boldsymbol{\theta}}={\boldsymbol{\theta}}_{0}$). We assume that $\mathcal{J}$ is twice continuously differentiable with respect to $({\boldsymbol{z}},{\boldsymbol{\theta}})$ and that the Hessian $\nabla_{{\boldsymbol{z}},{\boldsymbol{z}}}\mathcal{J}({\tilde{{\boldsymbol{z}}}},{\boldsymbol{\theta}}_{0})$ is positive definite. Applying the Implicit Function Theorem to (3.4) yields the existence of an operator ${\boldsymbol{z}}^{\star}:\mathcal{N}({\boldsymbol{\theta}}_{0})\to{\mathbb{R}}^{n_{z}}$, defined on a neighborhood of ${\boldsymbol{\theta}}_{0}$, such that

[TABLE]

for all ${\boldsymbol{\theta}}\in\mathcal{N}({\boldsymbol{\theta}}_{0})$. That is, the mapping ${\boldsymbol{z}}^{\star}({\boldsymbol{\theta}})$, from the model discrepancy parameters to the optimal decision variable, is uniquely defined in a neighborhood of the ROMCO solution. Furthermore, the Jacobian of ${\boldsymbol{z}}^{\star}({\boldsymbol{\theta}})$, referred to as the post-optimality sensitivity operator, is given by

[TABLE]

where ${\boldsymbol{H}}=\nabla_{{\boldsymbol{z}},{\boldsymbol{z}}}\mathcal{J}({\tilde{{\boldsymbol{z}}}},{\boldsymbol{\theta}}_{0})$ and ${\boldsymbol{B}}=\nabla_{{\boldsymbol{z}},{\boldsymbol{\theta}}}\mathcal{J}({\tilde{{\boldsymbol{z}}}},{\boldsymbol{\theta}}_{0})$ are the $({\boldsymbol{z}},{\boldsymbol{z}})$ and $({\boldsymbol{z}},{\boldsymbol{\theta}})$ Hessians of $\mathcal{J}$, evaluated at $({\tilde{{\boldsymbol{z}}}},{\boldsymbol{\theta}}_{0})$. The linear approximation

[TABLE]

provides a computationally efficient approach to update the ROMCO solution. Specifically, the decision variables can be updated using $N$ FOM evaluations by computing (3.5), where ${\boldsymbol{\theta}}=\overline{{\boldsymbol{\theta}}}$ is the posterior mean arising from the Bayesian calibration. Uncertainty in this update can be quantified by computing posterior samples for ${\boldsymbol{\theta}}$ and pushing them forward through (3.5) to compute posterior decision variable samples.

We summarize our HDSA-MD approach in Figure 1.

4 Atmospheric contaminant control using operator inference

This section presents a control problem for a two-species advection-diffusion-reaction system. Let $u_{1}({\boldsymbol{x}},t)$ denote a contaminant that is to be neutralized by deploying a decontaminant $u_{2}({\boldsymbol{x}},t)$ within a spatial domain $\Omega$. Both $u_{1}$ and $u_{2}$ diffuse homogeneously and advect along a stationary and spatially heterogeneous velocity field ${\boldsymbol{v}}$, which is uncertain (later parameterized by ${\boldsymbol{\xi}}$). The following system of PDEs describes this scenario:

[TABLE]

Here, $T=0.4$ is the final time, $\kappa_{1}=\kappa_{2}=0.1$ are diffusion coefficients, and $\mathring{u}_{1}$ is the initial contaminant profile. The homogeneous Neumann boundary conditions eq. 4.1c prescribe zero-flux boundaries. The decontaminant enters the domain through the source term $s$ and reacts with the contaminant with rate $\rho=2$.

We consider a two-dimensional domain $\Omega\subset(0,1.2)^{2}$ that excludes obstacles and discretize in space with finite elements using $n_{x}=10{,}702$ nodes. This leads to a system of ordinary differential equations

[TABLE]

where ${\boldsymbol{u}}_{1}(t),{\boldsymbol{u}}_{2}(t),{\boldsymbol{s}}(t),\mathring{{\boldsymbol{u}}}_{1}\in{\mathbb{R}}^{n_{x}}$ are, respectively, the spatial discretizations of $u_{1}$, $u_{2}$, $s$, and $\mathring{u}_{1}$; ${\boldsymbol{M}}\in{\mathbb{R}}^{n_{x}\times n_{x}}$ is the finite-element mass matrix. The matrices ${\boldsymbol{A}}_{1},{\boldsymbol{A}}_{2}\in{\mathbb{R}}^{n_{x}\times n_{x}}$ and $\mathbf{R}_{1},\mathbf{R}_{2}\in{\mathbb{R}}^{n_{x}\times n_{x}n_{x}}$ correspond to the discretization of the linear and bilinear operators in (4.1), respectively. The spatial domain, mesh, and initial state are shown in Figure 2.

4.1 Control problem formulation

The goal is to manipulate the source term ${\boldsymbol{s}}(t)$ to minimize the contaminant $u_{1}$ in a protection region that is represented by $\psi({\boldsymbol{x}})=\exp\left(-10\|{\boldsymbol{x}}-{\boldsymbol{x}}_{p}\|^{2}\right)$, where ${\boldsymbol{x}}_{p}=(0.9,0.35)$. We seek to minimize the contaminant concentration in the protection region, measured as $\|u_{1}({\boldsymbol{x}},t)\psi({\boldsymbol{x}})\|$, by injecting a decontaminant at $n_{q}=14$ discrete spatial points $\bar{{\boldsymbol{x}}}_{1},\ldots,\bar{{\boldsymbol{x}}}_{n_{q}}\in\Omega$. The protection region and injection locations are shown in Figure 2.

We model the source term as a sum of small-variance Gaussians centered at the injection points $\bar{{\boldsymbol{x}}}_{1},\ldots,\bar{{\boldsymbol{x}}}_{n_{q}}$,

[TABLE]

where $q_{i}(t;{\boldsymbol{z}})$, $i=1,2,\dots,n_{q}$, are functions of time parameterized by the decision variable ${\boldsymbol{z}}$. Note that the source magnitude is defined by $q_{i}(t;{\boldsymbol{z}})$ squared to ensure that the injection is nonnegative.

Letting ${\boldsymbol{u}}(t)=({\boldsymbol{u}}_{1}(t)^{\mathsf{T}}~{\boldsymbol{u}}_{2}(t)^{\mathsf{T}})^{\mathsf{T}}\in{\mathbb{R}}^{2n_{x}}$, $\boldsymbol{\psi}\in{\mathbb{R}}^{n_{x}}$, and ${\boldsymbol{q}}(t;{\boldsymbol{z}})=(q_{1}(t;{\boldsymbol{z}})^{\mathsf{T}},q_{2}(t;{\boldsymbol{z}})^{\mathsf{T}},\dots,q_{n_{q}}(t;{\boldsymbol{z}})^{\mathsf{T}})^{\mathsf{T}}\in{\mathbb{R}}^{n_{q}}$ denote the discretized state, protection zone weights, and injection magnitudes, respectively, we consider the optimization problem

[TABLE]

where $\odot$ denotes the Hadamard (elementwise) product and ${\boldsymbol{z}}\in{\mathbb{R}}^{n_{z}}$ defines the $n_{q}$ injection magnitudes at a collection of $n_{s}$ time nodes, $n_{z}=n_{q}n_{s}$. The regularization term $\gamma\|{\boldsymbol{q}}(t;{\boldsymbol{z}})\odot{\boldsymbol{q}}(t;{\boldsymbol{z}})\|^{2}=\gamma\sum_{i=1}^{n_{q}}q_{i}(t;{\boldsymbol{z}})^{4}$ induces a cost, scaled by the regularization coefficient $\gamma=10^{-5}$, to inject the decontaminant.333Note that the regularization involves the norm of ${\boldsymbol{q}}(t;{\boldsymbol{z}})\odot{\boldsymbol{q}}(t;{\boldsymbol{z}})$ since regularizing with the norm of ${\boldsymbol{q}}(t;{\boldsymbol{z}})$ would be equivalent to an optimization problem involving $\ell^{1}$ regularization and non-negativity constraints, thus leading to a lack of differentiability that hinders the post-optimality analysis.

Our goal is to compute (or approximate) the solution of (4.4) in real-time. However, the ODE system (4.2) depends on the velocity field ${\boldsymbol{v}}$, which is only known online when the optimization solution is required. For this experiment, the velocity field is parameterized by the vector ${\boldsymbol{\xi}}=(\xi_{1},\xi_{2},\xi_{3},\xi_{4})\in{\mathbb{R}}^{4}$ as

[TABLE]

We assume an offline training wind field parameter vector $\tilde{{\boldsymbol{\xi}}}=(1,1,10,75)^{\mathsf{T}}$. During the online phase, the testing wind field parameter vector is ${\boldsymbol{\xi}}^{\dagger}=(1.5,1.5,8,50)^{\mathsf{T}}$. The two wind fields are shown in Figure 3.

We will train an operator inference ROM offline using the training wind field parameter vector $\tilde{{\boldsymbol{\xi}}}$. In the online phase, we solve the optimal control problem using the ROM, compute a single full-model evaluation (4.2) with testing wind field parameter vector ${\boldsymbol{\xi}}^{\dagger}$, and update the optimal control using HDSA-MD. The subsections that follow detail these steps.

4.2 Operator inference ROM

Training data is generated by constructing $n_{c}=5$ smooth (as functions of time) random control vectors ${\boldsymbol{q}}^{(\ell)}(t)$, $\ell=1,\ldots,n_{c}$, and solving the full-order system eq. 4.2 with these controls and the training wind field parameter vector $\tilde{{\boldsymbol{\xi}}}$. The training controls are shown in the left panel of Figure 4. A proper orthogonal decomposition (POD) [39] is computed for the two states separately to produce basis matrices ${\boldsymbol{V}}_{1}\in{\mathbb{R}}^{n_{x}\times r_{1}}$ and ${\boldsymbol{V}}_{2}\in{\mathbb{R}}^{n_{x}\times r_{2}}$ so that the states are approximated as

[TABLE]

for $\widehat{{\boldsymbol{u}}}_{1}(t)\in{\mathbb{R}}^{r_{1}}$ and $\widehat{{\boldsymbol{u}}}_{2}(t)\in{\mathbb{R}}^{r_{2}}$. In our example, $r_{1}=19$ and $r_{2}=33$ are determined as the smallest basis dimensions such that the POD residual energy is less than $10^{-5}$; see (A.1) for details.

Inserting eq. 4.8 into eq. 4.2 and left multiplying the $k$th equation by ${\boldsymbol{V}}_{k}^{\mathsf{T}}$, we obtain a system of ODEs

[TABLE]

where the system matrices are given by

[TABLE]

Here, $\boldsymbol{\Phi}\in{\mathbb{R}}^{n_{x}\times n_{q}}$ corresponds to evaluating the injection basis functions $\phi_{i}$ from (4.3), $i=1,2,\dots,n_{q}$, at the spatial nodes.

Classical model-reduction methods use the matrices of the FOM eq. 4.2 to construct the ROM matrices eq. 4.10, but this requires intrusive access to the FOM code to explicitly construct ${\boldsymbol{M}},{\boldsymbol{A}}_{1},{\boldsymbol{A}}_{2},\mathbf{R}_{1},\mathbf{R}_{2},$ and $\boldsymbol{\Phi}$. We construct a ROM using operator inference [12, 21, 32], a non-intrusive alternative to classical methods. Specifically, operator inference postulates the model structure (4.9) and estimates the matrices (4.10) via least-squares regression using the training data snapshots. Details of the regression problem formulation and solution are given in Appendix A.

4.3 ROM-constrained optimization

Given our operator inference ROM, we consider the ROMCO problem,

[TABLE]

where $\widehat{{\boldsymbol{u}}}(t)=(\widehat{{\boldsymbol{u}}}_{1}(t)^{\mathsf{T}}~~\widehat{{\boldsymbol{u}}}_{2}(t)^{\mathsf{T}})^{\mathsf{T}}\in{\mathbb{R}}^{r_{1}+r_{2}}$.

The discretized control is defined at $n_{s}=100$ time nodes and $n_{q}=14$ spatial locations, hence the optimization algorithm must search a $n_{z}=1400$ dimensional space. The ROMCO problem eq. 4.11 is solved using derivative-based optimization, where objective-function gradient evaluations and Hessian-vector products are computed using the adjoint method. This enables efficient computation of derivatives in that the number of ROM evaluations per derivative evaluation is independent of the decision variable dimension $n_{z}$. The right panel of Figure 4 shows the optimized controls ${\boldsymbol{q}}(t;{\boldsymbol{z}})$. Observe that some components are nonzero over the entire time interval while others are only nonzero over a subinterval.

The controls shown in Figure 4 are optimal for the ROM constructed using the training wind field defined by ${\boldsymbol{\xi}}=\tilde{{\boldsymbol{\xi}}}$. We seek a solution in the test case when the wind field is defined by ${\boldsymbol{\xi}}={\boldsymbol{\xi}}^{\dagger}$. To assess the performance of the ROMCO problem solution, we simulate the system using the full-order model with the testing wind field and ROMCO controls. The state solution is depicted in Figure 5, where we observe that a nontrivial amount of contaminant reaches the protection zone.

4.4 Model discrepancy analysis

As seen in Figure 5, the ROMCO control solution is insufficient. We seek to improve the controls using only a single evaluation of the FOM solution operator ${\boldsymbol{S}}(t;\tilde{{\boldsymbol{z}}},{\boldsymbol{\xi}}^{\dagger})$. As outlined in Section 3.1, we use this single FOM evaluation to calibrate a representation of the discrepancy ${\boldsymbol{\delta}}(t;{\boldsymbol{z}},{\boldsymbol{\theta}})\approx{\boldsymbol{S}}(t;{\boldsymbol{z}},{\boldsymbol{\xi}}^{\dagger})-{\boldsymbol{V}}\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}})$, where ${\boldsymbol{V}}$ is the concatenation of the basis matrices ${\boldsymbol{V}}_{1}$ and ${\boldsymbol{V}}_{2}$, and $\widehat{{\boldsymbol{S}}}(t;{\boldsymbol{z}})$ is the solution operator for the ROM (4.9).

The decision-variable dimension is $n_{z}=1400$ and yet we have only $N=1$ FOM evaluation. Consequently, there is significant posterior uncertainty in the model discrepancy operator. However, notice that the objective function $\mathcal{J}({\boldsymbol{z}},{\boldsymbol{\theta}}_{0})$ has small curvature in many directions. Consider a Taylor approximation

[TABLE]

where the linear term is omitted since $\nabla_{\boldsymbol{z}}\mathcal{J}({\tilde{{\boldsymbol{z}}}},{\boldsymbol{\theta}}_{0})=\boldsymbol{0}$ and the higher-order (greater than quadratic) terms are neglected.

Let ${\boldsymbol{H}}\boldsymbol{w}_{i}=\lambda\boldsymbol{w}_{i}$ denote the eigenvalues and eigenvectors of the Hessian ${\boldsymbol{H}}$, $i=1,2,\dots,n_{z}$. We assume the eigenvalues are ordered as $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n_{z}}$. For a small eigenvalue $\lambda_{i}$, directions for which ${\boldsymbol{z}}-{\tilde{{\boldsymbol{z}}}}$ align with the eigenvector $\boldsymbol{w}_{i}$ will have a negligible $\mathcal{O}(\lambda_{i})$ magnitude influence on the objective function. However, observing (3.5), we see that ${\boldsymbol{\theta}}$’s such that ${\boldsymbol{B}}{\boldsymbol{\theta}}$ aligns with $\boldsymbol{w}_{i}$ will result in optimal solution updates that scale like $\lambda^{-1}$. Hence, if posterior samples align with eigenvectors corresponding to small eigenvalues of the Hessian, then we will observe large variances in optimal solution updates despite the fact that these large variances yield small changes in the objective function. To avoid such scenarios, we consider an alternative form of the optimal solution update that projects the post-optimality sensitivity onto the subspace of the Hessian’s leading eigenvectors.

Let ${\boldsymbol{P}}\in{\mathbb{R}}^{n_{z}\times n_{z}}$ be the orthogonal projector onto the span of the $r$ leading eigenvectors $\{\boldsymbol{w}_{j}\}_{j=1}^{r}$. The projected optimal solution update is given by

[TABLE]

In this example, we retain the $r=15$ leading eigenvectors of the Hessian as they exhibit two orders of magnitude decay, i.e., $\frac{\lambda_{15}}{\lambda_{1}}=\mathcal{O}(10^{-2})$.

The left panel of Figure 6 displays three components of the posterior mean controls alongside the optimized controls that were computed from the ROMCO problem eq. 4.11. These three components were chosen because they had the largest change (in the post-optimality update) amongst the $n_{q}=14$ components. The right panel of Figure 6 shows the time evolution of the protection region contaminant mass squared, $\|{\boldsymbol{u}}_{1}(t)\odot\boldsymbol{\psi}\|_{{\boldsymbol{M}}}^{2}$, which is the integrand in (4.4a). The reduction in contaminant demonstrates the benefits gained from using HDSA-MD with only a single full-order-model evaluation. Quantitatively, we can measure the control improvement by comparing the values of the FOMCO objective function. This requires an addition FOM evaluation that is typically not computed in practice, but done here for the purpose of validation. The FOMCO objective function, evaluated with the wind field parameter vector ${\boldsymbol{\xi}}^{\dagger}$ and controls defined by $\tilde{{\boldsymbol{z}}}$ and, is $0.355$. In comparison, the FOMCO objective function is $0.233$ when evaluated with the wind field parameter vector ${\boldsymbol{\xi}}^{\dagger}$ and controls defined by $\tilde{{\boldsymbol{z}}}-{\boldsymbol{P}}{\boldsymbol{H}}^{-1}{\boldsymbol{B}}\overline{{\boldsymbol{\theta}}}$, where $\overline{{\boldsymbol{\theta}}}$ denotes the posterior mean of the model discrepancy parameters. Hence, we realized a $34\%$ reduction in the FOMCO objective using only a single FOM evaluation to update the controls.

5 Wildfire ignition inversion using a neural operator ROM

This section presents an inverse problem for the ignition point of a wildfire in the Valles Caldera of northern New Mexico, which features heterogeneous topography and fuels. Additionally, the caldera is a high-consequence area due to nearby infrastructure, its status as a protected nature preserve, and a history of wildfire activity. Such inversion supports initialization of fire propagation codes as well as forensic analysis [38]. Wildfire modeling is complex as it involves a system of nonlinear PDEs that couples the atmosphere and fire propagation. We use the Fortran package WRF-SFIRE [25] to simulate the fire dynamics, as summarized below.

Let $u({\boldsymbol{x}},t)$ denote the level-set function, defined on a spatial domain $\Omega\subset{\mathbb{R}}^{2}$. Physical quantities resulting from the fire—such as temperature, heat flux, and the fire front—are computed from $u$ in a postprocessing phase. The interpretation of the level-set function is a signed distance to the fire front. Specifically, at a given time $t$, the points ${\boldsymbol{x}}\in\Omega$ where $u({\boldsymbol{x}},t)=0$ define the fire front, the points where $u({\boldsymbol{x}},t)<0$ define the burned region, and the points where $u({\boldsymbol{x}},t)>0$ define the unburned region. The magnitude of the level-set function measures distance from the fire front.

The wildfire is modeled using the coupled PDEs

[TABLE]

where $\boldsymbol{v}$ is three-dimensional atmospheric wind, $C>0$ is the model top of atmosphere, $\frac{D}{Dt}$ is the material derivative, and $T>0$ is the final time. The fire spread rate $\Psi$ depends on $\boldsymbol{v}$ as well as material properties such as fuel type, topographic slope, and fuel moisture. Equation eq. 5.1a is the level-set equation, while (5.1b) are the Euler equations for atmospheric wind $\boldsymbol{v}$. There is a two-way coupling between the fire and atmosphere. Winds will drive the fire in a particular direction via the dependence of $\Psi$ on $\boldsymbol{v}$, while the heat from the fire will create updraft that affects winds via the Euler equation forcing $f(u)$.

To find a solution to eq. 5.1, we must specify initial conditions and, as appropriate, boundary conditions. The initial condition for the level-set function $u$ is defined by the shifted pointwise Euclidean distance to the ignition point,

[TABLE]

where ${\boldsymbol{z}}\in{\mathbb{R}}^{2}$ is the fire ignition location and $c=-10$ is a shifting factor determined by the initial burn radius. Since $u$ is a signed distance function, no boundary conditions are imposed on the level-set function [25]. The initial condition for the wind velocity is user-specified based on knowledge of the atmosphere, and open boundary conditions are employed; see [40] for details.

The Fortran package WRF-SFIRE [25] handles all physical properties, empirical relationships, and numerics for eq. 5.1. The level-set function $u({\boldsymbol{x}},t)$ is discretized using a $173\times 133$ computational grid with spatial cells of size $60\times 60$ meters, resulting in $n_{u}=23009$ spatial nodes. We consider a time horizon of $T=8$ hours with a time-step of $2$ seconds. WRF-SFIRE employs a finite-volume spatial discretization with the $(x,y)$ mesh for eq. 5.1b being five times coarser than eq. 5.1a. Runge–Kutta time-steppers of order 2 and 3 are used for eq. 5.1a and eq. 5.1b, respectively. Operator splitting couples the equations, where $u$ is updated first and then provided as an input to $\boldsymbol{v}$.

5.1 Inverse problem formulation

Our goal is to estimate the ignition location ${\boldsymbol{z}}\in{\mathbb{R}}^{2}$ of the fire, given observations of the burn area at a limited number of time points $\{t_{j}^{\text{obs}}\}_{j=1}^{n_{\text{obs}}}$ that occur after ignition. We assume hourly observations that begin two hours after the fire ignition, i.e., observations at times $t_{j}^{\text{obs}}=1+j$ hours, $j=1,\dots,n_{\text{obs}}$, where $n_{\text{obs}}=7$.

In practice, satellite observations enable assessment of whether a spatial location has been burned, i.e., binary data indicating that a pixel is burned or unburned. Such data may be post-processed to estimate the fire front and compute the level-set function from it. Accordingly, we consider observational data of the form $\boldsymbol{Y}=({\boldsymbol{y}}_{1},\dots,{\boldsymbol{y}}_{n_{\text{obs}}})\in{\mathbb{R}}^{n_{u}\times n_{\text{obs}}}$, where $Y_{ij}$ is the observed signed-distance function at the $i^{th}$ spatial node ${\boldsymbol{x}}^{(i)}$ and time $t_{j}^{\text{obs}}$. Specifically, given the fire front $\Omega_{f}(t)\subset{\mathbb{R}}^{2}$ at time $t$, we denote the distance from a point $\boldsymbol{x}\in{\mathbb{R}}^{2}$ to the fire front as

[TABLE]

and define

[TABLE]

The ignition point is estimated by solving the optimization problem

[TABLE]

where ${\boldsymbol{u}}(t_{j}^{\text{obs}};{\boldsymbol{z}})\in{\mathbb{R}}^{n_{u}}$ denotes the solution of (5.1a) after spatial discretization, evaluated at time $t=t_{j}^{\text{obs}}$, with the level-set initial condition parameterized by ${\boldsymbol{z}}$ as defined in (5.2).

Solving (5.4) in real-time is intractable because the computational expense of WRF-SFIRE simulations limits the number of evaluations. Furthermore, it is impractical to implement sensitivity or adjoint equations due to the complex nonlinear feedbacks, operator splitting, and fire spread model characteristics used in WRF-SFIRE. Constructing an intrusive ROM would be challenging for the same reasons. Accordingly, we seek to evaluate WRF-SFIRE offline and train a non-intrusive ROM to approximate the level-set evolution.

5.2 Neural operator ROM

For the atmospheric transport problem in Section 4, we were able efficiently train a nonintrusive ROM via operator inference by leveraging the quadratic structure of the FOM. Because WRF-SFIRE involves more complex nonlinear structure and feedback mechanisms between the fire spread and atmospheric models, in this experiment we construct a nonintrusive ROM through neural operators.

Neural operators are a large class of methods that seek to leverage the strengths of neural networks to approximate mappings between infinite-dimensional spaces [20]. Many contemporary methods fall into this category, including deep operator networks (DeepONet) [23], Fourier neural operators (FNO) [22], flow map learning (FML) [10], and neural operators defined on POD bases [5]. These methods differ in the types of operators they approximate well. We employ FML because of our focus on predicting the level-set time evolution, along with the use of with neural operators defined on POD bases to facilitate learning from a limited number of training simulations. Recent work has extended FML from its original focus on autonomous dynamical systems [37] to include control inputs [35], parametric dependence [36], modal-space representations such as POD or Fourier series [42, 8, 9], and flow-map-constrained optimization [14].

In FML, it is frequently advantageous to learn the ROM on a coarse time grid since this can reduce the effect of error accumulation. Accordingly, we assume that the WRF-SFIRE simulation data is subsampled at $n_{\tau}=9$ points that correspond to the times $t=0,1,\dots,8$ hours. Letting $\mathring{{\boldsymbol{u}}}({\boldsymbol{z}})$ denote the spatial discretization of the initial condition parameterization (5.2), we consider the discrete-time dynamical system

[TABLE]

where ${\boldsymbol{f}}:{\mathbb{R}}^{n_{u}}\to{\mathbb{R}}^{n_{u}}$ corresponds to the integrating the WRF-SFIRE dynamics over a one-hour time interval. Note that ${\boldsymbol{f}}$ also depends on the atmospheric winds, although it is not explicit in the notation. Learning the coupling between the level-set function and the wind field from limited training simulations is challenging due to the complex nonlinear interactions between fire and the atmosphere. Instead, we will fix the atmospheric wind initial condition to a nominal estimate for all training data generation and only extract the level-set function as output from the WRF-SFIRE simulations. This corresponds, in our previous notation, to fixing ${\boldsymbol{\xi}}=\tilde{{\boldsymbol{\xi}}}$, where ${\boldsymbol{\xi}}$ represents the atmospheric wind initialization. Due to the fire spread model’s dependence on the transient wind field, (5.5) cannot represent the WRF-SFIRE dynamics exactly. However, with a fixed wind field initialization, (5.5) is a reasonable approximation to facilitate learning a ROM.

An initial wind field ${\boldsymbol{v}}({\boldsymbol{x}},0)=(0.1,0.1,0)^{\top}$ is fixed and the 20 simulations are run with different ignition locations to generate a training set (with 15 simulations) and validation set (with 5 simulations). The left panel of Figure 7 shows the terrain and ignition locations.

Our objective in FML is to approximate ${\boldsymbol{f}}$ in (5.5) using our training and validation set simulations. Since the dimension of a state time snapshot is $\mathcal{O}(10^{4})$, it is impractical to learn ${\boldsymbol{f}}:{\mathbb{R}}^{n_{u}}\to{\mathbb{R}}^{n_{u}}$ from only $20$ simulations. Rather, we project the dynamics into a POD basis and approximate the projected flow map. At least $r=100$ POD basis vectors are required to achieve a $1\%$ relative $\ell^{2}$ projection error on the validation trajectories.

Let ${\boldsymbol{V}}\in{\mathbb{R}}^{n_{u}\times r}$ denote the matrix containing the POD basis vectors. Multiplying eq. 5.5 by ${\boldsymbol{V}}^{\top}$ on the left motivates the FML-based reduced dynamical system

[TABLE]

where $\widehat{{\boldsymbol{f}}}({\boldsymbol{V}}^{\top}{\boldsymbol{u}}_{k})$ is a function we seek to learn that approximates the reduced WRF-SFIRE dynamics ${\boldsymbol{V}}^{\top}{\boldsymbol{f}}({\boldsymbol{V}}{\boldsymbol{V}}^{\top}{\boldsymbol{u}}_{k})$. Appendix B provides details regarding the network architecture of $\hat{{\boldsymbol{f}}}$, loss function, and training. The right panel of Figure 7 displays the fire front contour at $t=6$ hours for a validation set simulation and the corresponding ROM prediction.

5.3 ROM-constrained optimization

Once we have trained the ROM, the FOMCO problem eq. 5.4 is approximated using the ROMCO problem

[TABLE]

where $\widehat{{\boldsymbol{u}}}_{1+j}({\boldsymbol{z}})$, $j=1,2,\dots,n_{\text{obs}}$, satisfies eq. 5.6. To assess how well (5.7) recovers the ignition location, we generate a test set that consists of $5$ WRF-SFIRE simulations with ignition locations not used in the training and validation sets (Figure 7 depicts the ignition locations). Furthermore, the test set simulations use the initial wind field ${\boldsymbol{v}}({\boldsymbol{x}},0)=(4.33,2.5,0)^{\top}$, which has a significantly different magnitude and direction than the training and validation-set initial-wind fields. In the article’s notation, this test set wind field initialization corresponds to ${\boldsymbol{\xi}}^{\dagger}$.

For each testing dataset simulation, we consider the level-set from WRF-SFIRE as noise-free observations $\{{\boldsymbol{y}}_{j}\}_{j=1}^{n_{\text{obs}}}$ that enter the optimization objective. A trust-region method is used to solve (5.7) for each test set simulation. Gradients are computed by using automatic differentiation to compute the Jacobian of $\widehat{{\boldsymbol{f}}}$, and the Hessian is approximated using BFGS [29]. The data misfit—the difference between the test set data and the FOM solution using the estimated ignition location—is shown in the second column of Table 1. The rows correspond to each of the $5$ test cases. We observe misfit errors ranging from $0.1\%$ to $19\%$, with an average of $5.4\%$. The ignition location estimation error—the difference between the true and estimated ignition location—is shown in the second column of Table 2. The ignition location errors range from $281$ to $2145$ meters over the $5$ test cases, with an average of $967$ meters. This ignition location error is what we seek to reduce via HDSA-MD, which is described below and reported in the right columns of the tables.

5.4 Model discrepancy analysis

The large errors we observe in Table 1 and Table 2 highlight the deficiency of using a ROM trained with inaccurate atmospheric data. Nonetheless, each evaluation of WRF-SFIRE consists of taking over $10^{4}$ time steps, which requires considerable wall-clock-time, making it intractable to solve the FOMCO problem in real-time. Our goal is to get the best ignition estimate possible using only one WRF-SFIRE evaluation online.

For each of the 5 test cases, we generate a single FOM simulation using the test-set initial-wind field ${\boldsymbol{v}}({\boldsymbol{x}},0)=(4.33,2.5,0)^{\top}$ and ignition location given by the ROMCO solution ${\tilde{{\boldsymbol{z}}}}$ for the given test case. This single FOM evaluation is used to calibrate the model discrepancy. Let $\overline{{\boldsymbol{\theta}}}$ denote the posterior mean of the discrepancy calibration and $\overline{{\boldsymbol{z}}}=\tilde{{\boldsymbol{z}}}-{\boldsymbol{H}}^{-1}{\boldsymbol{B}}\overline{{\boldsymbol{\theta}}}$ denote the posterior mean ignition location. For each of the 5 test cases, we run WRF-SFIRE with the ignition location $\overline{{\boldsymbol{z}}}$ and test set initial wind field to measure the reduction in data misfit. The third columns of Table 1 and Table 2 display the posterior mean ignition location data misfit and ignition location estimation error, respectively, for each of the test cases. The fourth column of Table 1 and Table 2 shows the relative improvement from the ROMCO solution to the posterior mean ignition location. We observe $36\%-74\%$ reduction in the data misfit and $18\%-59\%$ reduction in the ignition location estimation error. The lowest percentage reduction corresponds to the test case where the error in $\tilde{{\boldsymbol{z}}}$ is already small. The consistent improvement attained using a single FOM evaluation in each test case demonstrates the benefit of HDSA-MD.

Lastly, we select one test case to demonstrate the utility of the model discrepancy analysis to provide uncertainty estimates. Figure 8 displays the posterior mean ignition location $\overline{{\boldsymbol{z}}}$ and a $95\%$ confidence ellipse around it (the posterior follows a Gaussian distribution), along with the ROMCO solution $\tilde{{\boldsymbol{z}}}$ and the “ground truth” ignition location ${\boldsymbol{z}}^{\dagger}$. We observe that the posterior mean ignition location lies approximately half-way between the ROMCO solution and the true ignition location, reflecting the approximately $50\%$ reduction in estimation error. Furthermore, the true ignition location ${\boldsymbol{z}}^{\dagger}$ lies within the $95\%$ confidence ellipse. Hence, the posterior distribution is well-calibrated despite the fact that it only utilized a single FOM evaluation.

6 Conclusion

Enabling optimization within wall-clock-time constraints is essential for developing digital twins that adapt to incoming data and inform decision-making. This article introduces a novel combination of non-intrusive ROMs and HDSA-MD to address the challenge of rapidly computing a reliable optimization solution. Our numerical results demonstrate ROM-constrained optimization for two different systems of nonlinear partial differential equations. In both cases, the ROMs are trained using inaccurate atmospheric data, which is subsequently corrected using HDSA-MD with a single FOM evaluation involving updated atmospheric data. This is representative of a range of digital-twin applications where model parameters vary with time, and hence the ROM is necessarily trained with inaccurate parameters.

In scenarios where it is feasible, retraining the ROM using the additional FOM evaluations and resolving the updated optimization problem may be superior. Our approach is closely related since HDSA-MD consists of (i) calibrating a model discrepancy representation to augment the ROM, and (ii) performing a post-optimality update to approximate the discrepancy augmented optimization solution. The advantage of HDSA-MD is that it avoids parameter tuning and additional computation required for ROM retraining, thus making the process automatic and real-time performant. Furthermore, HDSA-MD is effective with as little as one additional FOM evaluation. Lastly, HDSA-MD provides efficient estimates of uncertainty due to the model discrepancy.

Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research Field Work Proposal Number 23-02526, and the Laboratory Directed Research and Development program at Sandia National Laboratories. Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s National Nuclear Security Administration contract DE-NA0003525. This written work is authored by an employee of NTESS. The employee, not NTESS, owns the right, title and interest in and to the written work and is responsible for its contents. Any subjective views or opinions that might be expressed in the written work do not necessarily represent the views of the U.S. Government. The publisher acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this written work or allow others to do so, for U.S. Government purposes. The DOE will provide public access to results of federally sponsored research in accordance with the DOE Public Access Plan. SAND2025-10864O.

Appendices

​​ A Operator Inference ROM for Atmospheric Contaminant Modeling

Training data generation

We construct $n_{c}=5$ smooth random control vectors ${\boldsymbol{q}}^{(\ell)}(t)$, $\ell=1,\ldots,n_{c}$, by drawing values for ${\boldsymbol{q}}^{(\ell)}(\bar{t})$ at $\bar{t}=0$, $T/3$, $2T/3$, and $T$ from the uniform distribution over $[0,5]$ and interpolating these values with a monotone piecewise cubic spline [11]. The full-order system eq. 4.2 is then solved with ${\boldsymbol{q}}(t;{\boldsymbol{z}})={\boldsymbol{q}}^{(\ell)}(t)$, resulting in snapshot matrices ${\boldsymbol{U}}_{1}^{(\ell)},{\boldsymbol{U}}_{2}^{(\ell)}\in{\mathbb{R}}^{n_{x}\times n_{t}}$ for each $\ell=1,\ldots,n_{c}$. Here, the $j$-th column of ${\boldsymbol{U}}_{k}^{(\ell)}$ is the ${\boldsymbol{u}}_{k}(t_{j})$ (the $k$-th variable at $j$-th time) with forcing ${\boldsymbol{q}}^{(\ell)}$.

POD basis generation

The basis matrices ${\boldsymbol{V}}_{1},{\boldsymbol{V}}_{2}$ are obtained through weighted proper orthogonal decomposition (POD), which minimizes the reconstruction error of the snapshot data in the ${\boldsymbol{M}}$-norm [1, 4, 39]. Specifically, for $k=1,2$, we form the concatenated snapshot matrix ${\boldsymbol{U}}_{k}=({\boldsymbol{U}}_{k}^{(1)}~{\boldsymbol{U}}_{k}^{(2)}~\cdots~{\boldsymbol{U}}_{k}^{(n_{c})})\in{\mathbb{R}}^{n_{x}\times n_{c}n_{t}}$, compute the symmetric eigendecomposition ${\boldsymbol{U}}_{k}^{\mathsf{T}}{\boldsymbol{M}}{\boldsymbol{U}}_{k}=\boldsymbol{\Psi}_{k}\boldsymbol{\Sigma}_{k}^{2}\boldsymbol{\Psi}_{k}^{\mathsf{T}}$, and define ${\boldsymbol{V}}_{k}$ to be the first $r_{k}$ columns of ${\boldsymbol{U}}_{k}\boldsymbol{\Psi}_{k}\boldsymbol{\Sigma}_{k}^{-1}$. The reduced dimension $r_{k}$ is chosen based on the POD singular values—the diagonal values of $\boldsymbol{\Sigma}_{k}$, written $\sigma_{1}>\sigma_{2}>\cdots>\sigma_{n_{c}n_{t}}>0$—so that the so-called POD residual energy is less than $10^{-5}$, i.e.,

[TABLE]

In our experiment, the smallest $r_{1},r_{2}$ that satisfy this criterion are $r_{1}=19$ and $r_{2}=33$.

Operator inference regression problem

To generate the ROM, we formulate a linear regression problem to estimate the operators eq. 4.10. To this end, we project the snapshot matrices onto the bases ${\boldsymbol{V}}_{k}$, $k=1,2$, i.e., $(\widehat{{\boldsymbol{u}}}_{k}^{1}~\cdots~\widehat{{\boldsymbol{u}}}_{k}^{K})={\boldsymbol{V}}_{k}^{\mathsf{T}}{\boldsymbol{M}}{\boldsymbol{U}}_{k}$, with the number of columns $K=n_{c}n_{t}$ corresponding the product of the number of FOM evaluations $n_{c}$ and the number of time snapshots $n_{t}$. The operator inference regression problem requires access to time derivatives, which we denote as $\dot{\widehat{{\boldsymbol{u}}}}_{k}^{j}$. We approximate the time derivatives using sixth-order finite differences of the training states. Then we estimate the operators by independently solving the two linear regression problems

[TABLE]

To select appropriate regularization constants $\lambda_{1},\lambda_{2}\geq 0$, we take the approach of [26, 27, 33]: for each pair in a set of potential values $(\lambda_{1},\lambda_{2})$, we solve the operator inference regression problem eq. A.2 and evaluate the resulting ROM eq. 4.9 with the same initial conditions and controls as were used to generate the full-order data, then select the $(\lambda_{1},\lambda_{2})$ that minimizes the average reconstruction error. In this example, the process results in $\lambda_{1}=0.1$ and $\lambda_{2}=10$.

Remark on efficient objective function evaluation

The integrand in the objective eq. 4.11a can be evaluated efficiently by pre-computing norm weights for the reduced state. Let $\tilde{{\boldsymbol{M}}}\in{\mathbb{R}}^{n_{x}\times n_{x}}$ be the matrix with entries $\tilde{M}_{ij}=\psi_{i}M_{ij}\psi_{j}$, where the $M_{ij}$ are the entries of the mass matrix ${\boldsymbol{M}}$ and the $\psi_{i}$’s are the entries of $\boldsymbol{\psi}$. Then

[TABLE]

The matrix ${\boldsymbol{V}}_{1}^{\mathsf{T}}\tilde{{\boldsymbol{M}}}{\boldsymbol{V}}_{1}\in{\mathbb{R}}^{r_{1}\times r_{1}}$ can be computed once and stored for repeated use so that the cost of evaluating the objective scales with $r_{1}$ instead of $n_{x}$.

​​ B Neural operator ROM for Wildfire Modeling

We seek a neural network approximation $\widehat{{\boldsymbol{f}}}({\boldsymbol{V}}^{\top}{\boldsymbol{u}}\,;\boldsymbol{\gamma})\approx{\boldsymbol{V}}^{\top}{\boldsymbol{f}}({\boldsymbol{V}}{\boldsymbol{V}}^{\top}{\boldsymbol{u}})$, where $\boldsymbol{\gamma}$ are the network parameters (weights and biases). Let $\widehat{{\boldsymbol{u}}}$ denote the approximate reduced state coordinates generated by repeated composition of $\widehat{{\boldsymbol{f}}}$, as in eq. 5.6, where $\widehat{{\boldsymbol{u}}}$ and $\widehat{{\boldsymbol{f}}}$ depend parametrically on $\boldsymbol{\gamma}$.

Let $M=16$ denote the number of training dataset simulations and ${\boldsymbol{u}}^{(m)}_{k}$ denote the $m^{th}$ training dataset state solution at time $t_{k}$, $m=1,2,\dots,M$, $k=0,1,\dots,n_{\tau}$. To train $\widehat{{\boldsymbol{f}}}$, we seek to determine weights and biases $\boldsymbol{\gamma}$ that minimize the $P$-step recurrent loss

[TABLE]

where the innermost sum is assumed to be zero if $k+p>n_{\tau}$ and ${\boldsymbol{V}}$ is the POD basis matrix. Setting $1<P<n_{\tau}$ seeks to balance the long-time stability of $\widehat{{\boldsymbol{u}}}$ against the additional training costs of incorporating further compositions of $\widehat{{\boldsymbol{f}}}$ [14, 31].

The flow-map network has 4 hidden layers with 200 neurons per layer and a GELU activation function between layers; the output layer has an identity activation. Training is done with Tensorflow using the Adam optimizer. We train for 1000 epochs with an initial learning rate of $4\times 10^{-3}$ that decays smoothly to $10^{-3}$. At the end of training, the flow-map has a relative $\ell^{2}$ composition error of $2.3\%$ over the validation set.