=Paper=
{{Paper
|id=Vol-2964/article_202
|storemode=property
|title=Surrogate Modeling for Physical Systems with Preserved Properties and Adjustable Tradeoffs
|pdfUrl=https://ceur-ws.org/Vol-2964/article_202.pdf
|volume=Vol-2964
|authors=Randi Wang,Morad Behandish
|dblpUrl=https://dblp.org/rec/conf/aaaiss/WangB21
}}
==Surrogate Modeling for Physical Systems with Preserved Properties and Adjustable Tradeoffs==
https://ceur-ws.org/Vol-2964/article_202.pdf
Surrogate Modeling for Physical Systems
with Preserved Properties and Adjustable Tradeoffs
Randi Wang1 , Morad Behandish2 *
1
University of Wisconsin-Madison, 2 Palo Alto Research Center (PARC)
1
randi.wang@wisc.edu, 2 moradbeh@parc.com
Abstract
Determining the proper level of details to develop and solve
physical models is usually difficult when one encounters new
engineering problems. Such difficulty comes from how to
balance the time (simulation cost) and accuracy for the phys-
ical model simulation afterwards. We propose a framework
for automatic development of a family of surrogate models
of physical systems that provide flexible cost-accuracy trade-
offs to assist making such determinations. We present both a
model-based and a data-driven strategy to generate surrogate
models. The former starts from a high-fidelity model gener-
ated from first principles and applies a bottom-up model or- Figure 1: Displacement and thermal fields of a piston simu-
der reduction (MOR) that preserves stability and convergence lated by solving its finite element model
while providing a priori error bounds, although the resulting
reduced-order model may lose its interpretability. The latter
generates interpretable surrogate models by fitting artificial the geometry, etc. Resolving these details often leads to sig-
constitutive relations to a presupposed topological structure nificant computational time and cost that may turn out to be
using experimental or simulation data. For the latter, we use overkill to predict the desired properties up to an adequate
Tonti diagrams to systematically produce differential equa- accuracy. The tradeoffs are simply unknown until the model
tions from the assumed topological structure using algebraic
topological semantics that are common to various lumped-
is simulated many times under various assumptions and us-
parameter models (LPM). The parameter for the constitutive ing different LOGs.
relations are estimated using standard system identification High-fidelity models (HFMs) for DPM are commonly
algorithms. Our framework is compatible with various spa- specified by one or more partial differential equations
tial discretization schemes for distributed parameter models (PDEs) and initial/boundary conditions (ICs/BCs). To solve
(DPM), and can supports solving engineering problems in these equations, one often performs a discretization in space
different domains of physics. over a mesh, e.g., using finite difference [20], finite element
[3], and finite volume [12] or spectral methods [15], and se-
Introduction lects a time-stepping (e.g., Euler or Runge Kutta [8]) scheme
for integration. The mesh resolution or the number of ba-
Determining the appropriate level of granularity (LOG) to sis functions are determined by the user, to satisfy accuracy
model a physical system for computer-aided simulation is and stability requirements. The computational time and cost
a nontrivial problem, as cost-accuracy tradeoffs are rarely grows at a quadratic (in 2D) or cubic (in 3D) rate with de-
clear upfront. Generally, when encountered with a new prac- creasing mesh length, if not worse [6, 18]. For example, the
tical problem, physicists and engineers attempt to model the deformation and temperature fields of the piston in Fig. 1 is
system as accurately as possible, starting from the first prin- governed by two coupled PDEs describing the elastodynam-
ciples. For example, Fig. 1 illustrates a multiphysics simu- ics and transient heat transfer. A typical dynamic thermo-
lation of a piston operating in high pressure and tempera- elastic simulation for a system like this with geometric de-
ture conditions. Generally, to make such predictions accu- tails and heterogeneous materials may take anywhere from
rate, various details of the geometric and physical models minutes to hours or days, depending on the mesh resolution,
might be accounted for, such as the small geometric features integration time-step, and available computing power. Using
(e.g., holes and grooves), complex material distribution over such HFM and simulations for every part of a complex en-
* Corresponding Author gineering system (e.g., an entire automobile) is impractical
Copyright © 2021for this paper by its authors. Use permitted under and in most cases even unnecessary.
Creative Commons License Attribution 4.0 International (CC BY A common solution to this problem is to use reduced-
4.0) order models (ROM) that make compromises in accuracy
to achieve a better computational performance. However, it
is not easy to determine what kind of geometric or material
features may be omitted or simplified to retain the key prop-
erties one cares about. For example, decreasing the resolu-
tion across the model or “de-featuring” the mesh may result
in prohibitively large errors, while a much simpler model,
perhaps engineered with a careful consideration of what fea-
tures are important for the predictive outcome, may still cap-
ture the essential properties at a small fraction of the HFM
simulation cost. The domain insight or validation data for
such determinations are not always available. A systematic
and domain-agnostic approach that can be applied automat-
ically to generate flexible ROMs from a given HFM or data Figure 2: Framework for surrogate modeling
is missing. This paper aims to close that gap by providing
such a framework that further provides users with “knobs” to
adjust cost-accuracy tradeoffs, with guaranteed properties, are generally used when a “white box” HFM is available.
prior to committing to an LOG for simulation; to customize Bottom-up MOR techniques aim at finding a good approx-
interpretable structures that can be trained by data; and a imation of the HFM such that the approximation is numer-
combination of both, if applicable. ically efficient and stable while preserving certain physical
We present a framework (Fig. 2) that can systematically and numerical properties of the HFM [19, 4]. The MOR pro-
generate a family of ROMs for a system that span the cost- cess usually starts with a large system of ODEs and pro-
accuracy tradeoff spectrum. At the one end of this spectrum, duces a significantly smaller system of ODEs that is guaran-
one has the HFM described by a number of PDEs, semi- teed to predict a similar solution for given stimulus. Fig. 3
discretized (i.e., discretized in space but not in time) into (a) illustrates the key ideas in MOR. By contrast, the data-
a large system of coupled ordinary differential equations driven methods, including machine learning (ML) or regres-
(ODEs) or differential algebraic equations (DAEs). One may sion approaches to parameter estimation for a small system
also have access to data sets representing the behavior, i.e., of ODEs, are generally used if the HFM is a “black box”,
response of the HFM (or actual physical system) to certain where some input and output data sets are given, but the
stimuli such as ICs/BCs or source/sink terms. model’s inner working mechanism is not available.
The generated family of surrogate models may pro- The commonly-used MOR methods for linear time-
vide one or a combination of properties: (1) flexible cost- invariant (LTI) systems are the balanced truncation (BT)
accuracy tradeoffs (2) a priori error bounds (3) a priori guar- method [9, 26, 31] and the rational Krylov subspace (RKS)
antees of preserving properties that are critical for physical method [1, 2, 14, 16], which are based on projection. The
fidelity (4) interpretability. Equipped with this framework, principle of the BT method is to remove the states that have
physicists and engineers can quickly choose the LOG at weak controllability and observability [9]. The BT method
which the ROM can model the system and predict its be- is not time-efficient for large-scale models that have more
haviors most efficiently up to the required accuracy. than a few thousands of variables because finding out such
states is a time-consuming process. By contrast, the RKS
Contributions method matches several most significant terms of the Tayler
This paper has three major contributions: series expansion of the ROM transfer function, expanded
around carefully selected frequencies, to those of the HFM.
1. A MOR approach from the system and control theory Although RKS is more scalable than BT, it has several draw-
is adapted in developing ROMs for DPMs, using semi- backs; for example, it does not guarantee preservation of the
discretization (i.e., discretization in space, but not time) HFM’s stability (i.e., even if the HFM is stable, the ROM
to covert PDEs to a large system of ODEs/DAEs. may not be) and the Taylor series expansion frequencies
2. A domain-agnostic mechanism is proposed to automat- have to be selected manually, which is nontrivial.
ically translate a presupposed topological structure for To overcome the above limitations of RKS, researchers
LPMs to a system of ODEs, which supports both single have developed an improved method based on an iterative
physics and multiphysics couplings. algorithm [17], which allows adaptively choosing the ex-
3. A hybrid approach for surrogate modeling is proposed, pansion frequencies using the mirror image of the poles of
which, not only maintains all desirable properties of the the obtained transfer function of ROM. This algorithm is
model-based method, but also retrieves the interpretability considered as a “gold standard” among the projection-based
of the physical system through imposing LPM structure. MOR methods that minimize the norm of approximation er-
ror between ROM and HFM transfer functions for a given
target model order1 . Compared to RKS, the IRKA has the
Related work advantage of the ability to automatically choose expansion
Different approaches to surrogate modeling of physical sys-
tems can be broadly categorized into model-based and data- 1
The “order” of a model is the number of state variables, and
driven methods. The model-based methods such as MOR the number of first-order ODEs, in the state-space representation.
frequencies by iteration while maintaining the model sta- Methods and Algorithms
bility [17]. The algorithm of IRKA is summarized in Fig. Below, we introduce our model-based and data-driven ap-
3 (b). However, it can neither guarantee monotonically de- proaches to bottom-up MOR of white box systems and top-
creasing errors with every iteration, nor ensure that the error down surrogate modeling of black box systems, respectively.
converges to a local minimum [5]. The CUmulative REduc- At the end, we provide guidelines on how to use both ap-
tion (CURE) scheme [22] was later developed to improve proaches in a hybrid setting for “gray box” systems.
several critical properties of the IRKA. The algorithm of
the CURE scheme can be found in [22] and is repeated in A Model-Based Approach
Fig. 3 (c). Particularly, to maintain the model stability, a For the model-based approach, we use the SPARK+CURE
stability-preserving, adaptive rational Krylov (SPARK) al- [22] due to its rigorous mathematical underpinnings that
gorithm was developed [22], which is usually embedded lead to strong guarantees for numerical simulation. In
in the CURE scheme to generate a family of stable ROMs a nutshell, this method has several important properties,
with increasing model orders by sequential accumulation in namely: (1) automatic discovery of proper expansion fre-
a single MOR process. Both IRKA and SPARK share the quencies; (2) guaranteed preservation of stability; (3) guar-
same model-reduction principle as the RKS methods and are anteed H2 −error convergence; and (4) an a priori H2 −error
both numerically efficient. However, compared to the IRKA bound.2 More specifically, the method not only allows adap-
method, the SPARK+CURE has further advantages such as tively choosing the expansion frequencies using the poles of
ensuring stability, monotonic convergence, and a priori er- the obtained ROM (e.g., similarly to IRKA), but also enables
ror guarantees, as well as automatic model order decision incrementally increasing the model order by cumulatively
making capabilities. updating the ROM error transfer function. Particularly, the
Table 1 illustrates a comparison of the properties of the error of the obtained ROM monotonically converges to zero
four MOR methods mentioned above. The major drawback (i.e., ROM converges to HFM) in the accumulation process.
of all of these projection-based MOR techniques is the loss The SPARK+CURE scheme can also receive a tolerance for
of physics-based structure and interpretability. There ex- error, based on which a termination criterion for cumulative
ist structure-preserving (e.g., symplectic) MOR approaches reduction can be defined. In other words, the user provides
[24] that can remedy this shortcoming; however, the added not only the HFM matrices but also a maximal value of er-
interpretability may come at the expense of losing the de- ror that can be overlooked, and the algorithm generates the
sirable properties in Table 1. We will present a simple hy- lowest-order ROM (hence the fastest to simulate) that is still
bridization strategy to augment CURE with a data-driven guaranteed to remain within the error tolerance, without ac-
approach to retrieve the physical structure with a small com- tually performing any numerical simulation.
promise in the a priori error guarantees.
A Data-Driven Approach
Table 1: Comparison of properties of four commonly-used
MOR methods (CURE is uniquely positioned) In a recent article [34], common reference semantics for
lumped parameter system modeling were presented, based
on algebraic topological foundations of network theory
[27, 7], which can serve as a unifying abstraction of system
modeling languages such as Modelica [11], Simulink [10],
linear graphs [28], and bond graphs [23], etc. A key advan-
tage of using this abstraction is the ability to automatically
map a given topological structure for the LPM (e.g., a cir-
cuit graph or mass/spring/damper network) to a set of gov-
erning ODEs. These ODEs have built-in conservation laws
in the LPM context such as Kirchhoff’s current and volt-
age laws for electrical and thermal circuits, superposition of
forces and Newton’s laws of motion in multi-body dynam-
ics, and so on. They also include constitutive laws associ-
Popular data-driven methods include, but are not lim- ated with lumped components such as springs and dampers
ited to, symbolic regression [29, 30], recurrent neural net- in mechanical systems, resistors and capacitors in analog cir-
works [32], evolved regulatory networks [13], and physics- cuits, conductors in heat transfer, and so on. The recipe for
informed neural networks [25]. These methods have demon- generating the ODEs from system topology in [34] is given
strated the ability to accurately replicate HFM after suffi- by Tonti diagrams [33] of network theory (Fig. 4 (c)). The
cient training and testing. Their development often requires Tonti diagram is a composition of topological and algebraic
domain-specific insight to select the proper set of differen- operators that map data associated with different cells in an
tial equations and directly build these equations into a con- oriented cell complex representation of the LPM network;
strained learning structure and/or penalize the loss function for instance, in an electrical circuit, the superposition of in-
by residual errors. In this paper, however, we provide an ap- coming/outgoing currents on incident wires (i.e., 1−cells) to
proach that only requires the user insight in choosing an ap-
propriate LPM topology while the system equations can be 2
The H2 −error is defined as the L2 −norm of the error transfer
automatically generated in a domain-agnostic fashion. function over the imaginary line in the frequency domain.
Figure 3: The workflow for the model-based approach to LPM construction. The user provides (a) the HFM, obtained by
semi-discretization of a set of PDEs into a large system of ODEs/DAEs (e.g., via finite element spatial discretization). The
algorithm re-arranges the second-order equations into twice as many first-order equations in the state-space representation. The
state-space representation can be projected to a lower-dimensional space by various MOR techniques, e.g., (b) IRKA. We use
an improved version of IRKA, called (c) the CURE scheme. The diagram in (c) is adopted from [22] with modification.
a junction (i.e., 0−cells) is captured by a boundary operator use ordinary least squares regression—although other ob-
(from 1−cells to 0−cells), whereas resistance, capacitance, jective functions and optimization techniques are certainly
inductance, and other constitutive relations are in-place al- applicable—to iteratively update the constitutive parameters
gebraic relations that keep data on 1−cells. These operators until the solution of the state equation fits the given data.
are represented by different types of arrows on the Tonti di- In a later section, we will apply this method to obtain
agram. The key advatange of using this approach to equa- an LPM for a mechanical problem (single physics) and a
tion generation is its generalizability to various domains of thermo-mechanical problem (coupled multiphysics).
physics and possible multi-physics, as the underlying topo-
logical and algebraic operations are common to mechanical, A Hybrid Approach
electrical, thermal, and other systems [33]. Both of the model-based and data-driven approached men-
Constitutive relationships in LPM capture “effective” tioned above have pros and cons; in particular:
phenomenological properties of the system at a certain LOG • The model-based approach has the advantage of starting
of choice, unlike the constitutive relations in DPM de- from first principles and requiring no data other than basic
rived directly from well-documented material properties. material properties that are used in the HFM (e.g., consti-
Except in cases where an LPM is directly generated from tutive part of the PDEs).
a system model with modular components (e.g., actual
• The model-based method also provides rigorous guaran-
springs/dampers in an automobile suspension assembly), the
tees that are rarely avialble in data-driven methods; how-
parameters for the artificial LPM components are not eas-
ever, the resulting ROMs are often not interpretable.
ily obtained from geometric and material properties. These
parameters must be estimated from data by solving an opti- • The data-driven method provides a mechanism to specify
mization problem (e.g., least squares regression). the desired LPM with interpretable constitutive relations
Figure 4 illustrates the workflow for the data-driven LPM (“artificial” components as 1−cells in the cell complex) if
construction. Given an experimental or simulation data set the user has insight to choose an appropriate structure.
and an LPM topology in Fig. 4 (a), we first convert the In particular, even though the SPARK+CURE method
LPM from the domain-specific format (e.g., Modelica) to has several useful properties, it has the following limita-
the domain-agnostic canonical form (i.e., oriented cell com- tions: (1) the HFM itself must be dissipative (hence stable)
plex) as shown in Fig. 4 (b), using the semantics provided in to begin with (2) the resulting HFM is a system of first-
[34]. Each 1−cell is associated with a symbolic constitutive order ODEs/DAEs with dense state-space matrices which
relation and given an initial value to the constitutive parame- may not be refactored into a second-order system (e.g., with
ter. After selecting a state variable of interest, the state equa- mass, spring, and damper matrices in mechanical LPM) for
tions (i.e., system of second-order ODEs/DAEs) are gen- component-wise interpretability and (3) the rigorous guar-
erated by tracing groups of paths along the Tonti diagram antees are only valid for LTI systems, although there are
of network theory with the appropriate physical types [34]. several ways in which it can be generalized to nonlinear sys-
There are a total of 8 different options for state variables, tems, some of which are ongoing research.
each of which corresponds to a different groups of paths On the other hand, the regression-based system identifi-
[34]. Once the system of ODEs/DAEs are assembled, we cation for data-driven fitting becomes impractical when the
Figure 4: The workflow for the data-driven approach to LPM construction. The user provides (a) a topology for the LPM (i.e.,
symbolic network of inter-connected components), which is then converted to (b) the common language of abstract oriented
cell complexes [34]. (c) The Tonti diagram [33] converts the cell complex representation to (d) a system of symbolic ODEs
with unknown constitutive parameters. The parameters are learned from data using standard system identification techniques.
HFM is too costly to simulate (to obtain synthetic data) and other hand, comes into play for numerical simulation of the
the experimental measurements of adequate quality are not system of ODEs/DAEs (before or after MOR) by finite time-
available. One can always use as many data points as one stepping to approximate integration with given ICs.
can obtain within the computational and experimental bud-
get, but the lack of guarantees in predicting out-of-training Mechanical and Thermal Examples
inputs undermines the ROM’s reliability.
We will generate data using two geometric domains; namely,
To remedy the shortcomings of the the model-based and
a simple cylinder (Fig. 6) to illustrate the basic ideas with a
data-driven approaches, we developed a “hyrbid” approach
closed-form ground truth to compare against, and a piston
to achieve the best of both worlds. For each generated ROM
assembly (Fig. 8) to demonstrate a slightly more compli-
M from the SPARK+CURE method, we use this ROM to
cated case involving one-way coupled thermo-elasticity. In
produce training data for another surrogate LPM M0 of the
both examples, we assume linear elastic material properties
same order r n via system identification. Although the
undergoing dynamic mechanical and/or thermal loads. The
training data is itself erroneous, its H2 −error is bounded by
BCs in these examples are fixed displacement at some sur-
the CURE framework (denoted by εM ). The error between
faces and uniform pressure and heat fluxes at other surfaces.
the two ROMs (denoted by ε̄rel ), on the other hand, can be
In principle, the approach can be applied to domains of ar-
computed by solving two small-scale algebraic Lyapunov
bitrary shapes, material distributions, ICs/BCs, and excita-
equations.3 Hence, the H2 −error for the surrogate LPM can
tions. More extensive testing with strongly (e.g., two-way)
also be guaranteed via a triangular inequality:
coupled mutliphysics problems is required to further vali-
ε̄M0 ≤ ε̄M + ε̄rel . (1) date the approach, which is suitable for a full paper.
Applications and Results Preliminary Results
Spatial Discretization Consider the homogeneous cylinder in Fig. 6 (a) with one
To generate data for our data-driven method by using the end fixed to the ground and the other end bearing a constant
ROM generated from the SPARK+CURE method, we will unit pressure. The cylinder is discretized using second-order
use LTI models and FEA discretization [3], but virtually tetrahedral finite elements, leading to 17,430 nodal displace-
every spatial discretization works as long as it generates a ment variables. Using the vertical average displacement of
system of LTI ODEs/DAEs. Note that we do not perform the top surface of the cylinder as the solution of interest, we
temporal discretization for MOR, as the SPARK+CURE use the SPARK+CURE to generate a family of ROMs. Fig-
method is best described in terms of continuous and differ- ure 5 (a) shows the a priori relative H2 −error bounds for
entiable functions in time. Temporal discretization, on the each of them. It can be observed that the error bound mono-
tonically decreases with increasing order of the ROMs. For
3 visualization, we compare the simulation results of the ROM
Note that solving such a system for the full-order HFM would
be almost as prohibitive as numerical simulation. and the original HFM (finite elements simulation) as shown
(a) Cylinder (b) Mesh (c) Cell complex (d) mass-spring-
damper system
(a) A priori relative H2 −error bound
(e) Simulation results comparison
Figure 6: The topological structure of the physical space of
a lumped mass-spring-damper system and the comparison
of simulation results between the HFM, the reduced FEA
model, and the optimal model
(b) Simulation results comparison
(NRMSE)4 of 4.52%. The relative H2 −error bound com-
Figure 5: Simulation results comparison between HFM and puted from (1) against HFM is 0.184.
a family of ROMs, and a priori relative H2 −error bounds. Next, we apply the approach to a slightly more complex
problem with weak thermo-elastic coupling, over a piston
geometry shown in Fig. 8 (a), in which the temperature
change has an impact on structural field, not vice versa. Parts
of the piston are combined into one single part to avoid rela-
in Fig. 5 (b). We select the ROM of order r = 50 to generate tive motion. The piston bears a unit pressure (red arrows) on
simulation data, whose a priori relative H2 −error bound is the top surface, a constant heat flux (blue pins) on the pis-
10−2.4509 ≈ 0.35% as shown in Fig. 5 (a). ton head, and the crank is fixed (green pins). The solution of
interest is the vertical average displacement of the top sur-
To find an interpretable ROM with a desired topologi- face. The data is obtained by simulating a ROM generated
cal structure, we initialize an oriented cell complex (Fig. from the SPARK+CURE method, with an a priori relative
6 (c)) representing a lumped mass-spring-damper network H2 −error bound of 0.0031, computed from (1).
with 2 degrees of freedom (Fig. 6 (d)). We label the 1−cells To retrieve interpretability, we initialize a pair of con-
with L3 , L6 for masses, L2 , L5 for springs, and L1 , L4 for nected oriented cell complexes representing a mechanical
dampers. The Tonti diagram of mechanical LPM is shown in LPM (mass-spring-damper network) connected to a thermal
Fig. 7a, and the paths traced to generate the system of ODEs LPM (resistance-conductance network) by a transformer
are shown in Fig. 7 (b). We apply an ordinary least squares (TF) (Fig. 11), where we label the 1−cells with L3 , L6
regression [21] to find the values of the lumped masses, stiff- for masses, L2 , L5 for springs, L1 , L4 for dampers, L7 ,
ness coefficients, and damping coefficients and compare the L10 , L11 for thermal resistors, and L8 , L9 for thermal con-
simulation result for the optimal LPM against the ROM and ductors. The Tonti diagram of generalized network systems
the HFM (Fig. 6 (e)). It can be observed that the simula-
tion results between the ROM and the HRM are close and 4
The NRMSE is defined as the L2 −norm of the error signal
the simulation results between the optimal LPM and the (between ROM and LPM) in the time domain, normalized by the
ROM match well, with a normalized root-mean-square error variation interval of the ROM.
(a) Tonti diagram (b) Paths
Figure 7: Tonti diagram of lumped mass-spring-damper sys- Figure 10: A thermo-mechanical LPM with a TF.
tem (with lumped sources) and paths for generating govern-
ing equations from LPM topology.
Figure 11: The topological structure of the physical space
of a thermo-mechanical LPM with a transformer (TF) coou-
pling the two physics.
(a) Piston (b) Mesh of piston
Figure 8: A Piston and its FEA mesh for HFM simulation. Conclusions
We proposed a framework for automatically generating a
family of surrogate models of physical systems with dif-
ferent cost-accuracy tradeoffs. We presented model-based,
data-driven, and hybrid techniques that provide a priori guar-
antees of preserved properties including stability, conver-
gence, and error bounds, as well as an ability to enforce an
interpretable topological structure while assessing its impact
on the error bound at a small computational cost. In particu-
lar, we used the SPARK+CURE algorithm to develop a first
instance of a ROM for rapid simulation and data generation.
We use this ROM to train another surrogate ROM (inter-
pretable LPM) and measure the additional errors by solving
algebraic Lyapunov equations. To systematically develop an
(a) Generalized Tonti diagram (b) Paths
interpretable LPM, we used Tonti diagrams of network the-
ory to generate systems of ODEs from the presupposed topo-
Figure 9: Generalized Tonti diagram of multi-domain LPM logical structure and system identification to estimate the un-
and paths for generating governing equations. known constitutive parameters. This approach avoids simu-
lating the HFM and enables quantifying how well the as-
sumed topology fits the ROM.
[34] is shown in Fig. 9 (a) and the paths traced to generate This framework can assist researchers and engineers in
ODEs of the multi-domain lumped parameter systems are making decisions on the proper LOG to develop physical
shown in Fig. 9 (b). We apply an ordinary least squares re- models of new engineering problems, and will open up new
gression to obtain the optimal values of the lumped mass, research directions. Possible areas of further research in-
stiffness coefficients, damping coefficients, thermal conduc- clude extension to nonlinear systems using methods such as
tance, and thermal resistance. We compare the solutions piecewise linearization, dynamic mode decomposition, and
of the HFM, the ROM and the optimal LPM in Fig. 12, Koopman operator, and applying the ROM models to solve
where the NRMSE of the parameter estimation is 3.11% inverse problems (e.g., engineering design).
and the relative H2 −error bound computed from (1) against
the HFM is 0.0124. Note that the mechanical vibration time Acknowledgments
scale in this case is much faster than the heat diffusion, so the
displacement reaches steady-state value of −1.393 × 10−4 This material is based upon work supported by the De-
meters early on. The displacement visible in the figure is fense Advanced Research Projects Agency (DARPA) under
caused by thermal expansion. Agreements No. HR00111990029 and HR00112090065.
[11] Elmqvist, H.; and Mattsson, S. E. 1997. An introduc-
tion to the physical modeling language Modelica. In
Proceedings of the 9th European Simulation Sympo-
sium, ESS, volume 97, 19–23. Citeseer.
[12] Eymard, R.; Gallouët, T.; and Herbin, R. 2000. Finite
volume methods. Handbook of numerical analysis 7:
713–1018.
[13] François, P.; Hakim, V.; and Siggia, E. D. 2007. Deriv-
ing structure from evolution: metazoan segmentation.
Molecular systems biology 3(1).
[14] Freund, R. W. 2004. SPRIM: structure-preserving
reduced-order interconnect macromodeling. In
IEEE/ACM International Conference on Computer
Aided Design, 2004. ICCAD-2004., 80–87. IEEE.
Figure 12: Comparsion of simulation results between the re- [15] Gottlieb, D.; and Orszag, S. A. 1977. Numerical
duced FEA model and the algebraic topological model of a analysis of spectral methods: theory and applications.
thermo-mechanical system. SIAM.
[16] Grimme, E. 1997. Krylov projection methods for
References model reduction. Ph.D. thesis.
[1] Bai, Z. 2002. Krylov subspace techniques for reduced- [17] Gugercin, S.; Antoulas, A. C.; and Beattie, C. 2008.
order modeling of large-scale dynamical systems. Ap- H2 model reduction for large-scale linear dynamical
plied numerical mathematics 43(1-2): 9–44. systems. SIAM journal on matrix analysis and appli-
[2] Bai, Z.; and Su, Y. 2005. Dimension reduction of large- cations 30(2): 609–638.
scale second-order dynamical systems via a second- [18] Inc., L. 2020. Understanding Mesh Refinement and
order Arnoldi method. SIAM Journal on Scientific Conformal Mesh in FDTD. https://support.lumerical.
Computing 26(5): 1692–1709. com/hc/en-us/articles/360034382594-Understanding-
[3] Bathe, K.-J. 2006. Finite element procedures. Klaus- Mesh-Refinement-and-Conformal-Mesh-in-FDTD.
Jurgen Bathe. [19] Lohmann, B.; and Salimbahrami, B. 2000. Introduc-
[4] Baur, U.; Benner, P.; and Feng, L. 2014. Model tion to Krylov subspace methods in model order reduc-
order reduction for linear and nonlinear systems: a tion. Methods and applications in automation 1–13.
system-theoretic perspective. Archives of Computa- [20] Mazumder, S. 2015. Numerical methods for partial
tional Methods in Engineering 21(4): 331–358. differential equations: finite difference and finite vol-
[5] Beattie, C. A.; and Gugercin, S. 2009. A trust region ume methods. Academic Press.
method for optimal h 2 model reduction. In Proceed-
[21] Miller, S. J. 2006. The method of least squares. Math-
ings of the 48h IEEE Conference on Decision and Con-
ematics Department Brown University 8: 1–7.
trol (CDC) held jointly with 2009 28th Chinese Control
Conference, 5370–5375. IEEE. [22] Panzer, H. K. 2014. Model order reduction by Krylov
[6] Boz, Z.; Erdogdu, F.; and Tutar, M. 2014. Effects of subspace methods with global error bounds and auto-
mesh refinement, time step size and numerical scheme matic choice of parameters. Ph.D. thesis, Technische
on the computational modeling of temperature evo- Universität München.
lution during natural-convection heating. Journal of [23] Paynter, H. M. 1961. Analysis and design of engineer-
Food Engineering 123: 8–16. ing systems. MIT press.
[7] Branin, F. H. 1966. The algebraic-topological basis for [24] Peng, L.; and Mohseni, K. 2016. Symplectic model
network analogies and the vector calculus. In Sympo- reduction of Hamiltonian systems. SIAM Journal on
sium on generalized networks, 453–491. Scientific Computing 38(1): A1–A27.
[8] Butcher, J. C. 1987. The numerical analysis of ordi- [25] Raissi, M.; Perdikaris, P.; and Karniadakis, G. E. 2019.
nary differential equations: Runge-Kutta and general Physics-informed neural networks: A deep learning
linear methods. Wiley-Interscience. framework for solving forward and inverse problems
[9] Chahlaoui, Y.; Gallivan, K. A.; Vandendorpe, A.; and involving nonlinear partial differential equations. Jour-
Van Dooren, P. 2005. Model reduction of second-order nal of Computational Physics 378: 686–707.
systems. In Dimension Reduction of Large-Scale Sys- [26] Reis, T.; and Stykel, T. 2008. A survey on model re-
tems, 149–172. Springer. duction of coupled systems. In Model order reduc-
[10] Chaturvedi, D. K. 2009. Modeling and simulation of tion: theory, research aspects and applications, 133–
systems using MATLAB and Simulink. CRC Press. 155. Springer.
[27] Roth, J. P. 1955. An application of algebraic topology
to numerical analysis: On the existence of a solution
to the network problem. Proceedings of the National
Academy of Sciences 41(7): 518–521.
[28] Rowell, D.; and Wormley, D. N. 1997. System dynam-
ics: an introduction. Prentice Hall.
[29] Schmidt, M.; and Lipson, H. 2009. Distilling free-
form natural laws from experimental data. science
324(5923): 81–85.
[30] Schmidt, M. D.; Vallabhajosyula, R. R.; Jenkins, J. W.;
Hood, J. E.; Soni, A. S.; Wikswo, J. P.; and Lipson, H.
2011. Automated refinement and inference of analyt-
ical models for metabolic networks. Physical biology
8(5): 055011.
[31] Stykel, T. 2002. Analysis and numerical solution of
generalized Lyapunov equations. Institut für Mathe-
matik, Technische Universität, Berlin .
[32] Sussillo, D.; and Abbott, L. F. 2009. Generating coher-
ent patterns of activity from chaotic neural networks.
Neuron 63(4): 544–557.
[33] Tonti, E. 2013. The mathematical structure of classical
and relativistic physics. Springer.
[34] Wang, R.; and Shapiro, V. 2019. Topological semantics
for lumped parameter systems modeling. Advanced
Engineering Informatics 42: 100958.