## Abstract

Engineering design often involves representation in at least two levels of abstraction: the system-level, represented by lumped parameter models (LPMs), and the geometric-level, represented by distributed parameter models (DPMs). Functional design innovation commonly occurs at the system-level, followed by a geometric-level realization of functional LPM components. However, comparing these two levels in terms of behavioral outcomes can be challenging and time-consuming, leading to delays in design translations between system and mechanical engineers. In this paper, we propose a simulation-free scheme that compares LPMs and spatially discretized DPMs based on their model specifications and behaviors of interest, regardless of modeling languages and numerical methods. We adopt a model order reduction (MOR) technique that a priori guarantees accuracy, stability, and convergence to improve the computational efficiency of large-scale models. Our approach is demonstrated through the model consistency analysis of several mechanical designs, showing its validity, efficiency, and generality. Our method provides a systematic way to compare system-level and geometric-level designs, improving reliability and facilitating design translation.

## 1 Introduction

In the field of engineering design, the initial step toward creating innovative designs often involves system design, which entails conceptualizing and developing functional structures. Once a system is designed, the subsequent step involves developing its 3D geometry for manufacturing purposes. The objective of the system-based geometric design process is to identify at least one geometric realization that matches the system design. In this paper, we restrict our attention to mechanical systems whose representation at the system-level can be given by mass–spring–damper networks.

### 1.1 Motivation.

The distributed parameter models (DPMs; e.g., 3D geometric assemblies) and lumped parameter models (LPMs; e.g., mechanical mass–spring–damper networks) that describe the same physical system at different levels of abstraction are represented using different languages and semantics that cannot be directly translated into one another. For instance, Dassault Systémes offers two commercial software, dymola for system modeling and solidworks for geometric modeling. The models created in these programs are incompatible and cannot be automatically translated into each other owing to disparities in model types and representations [1–4]. This gap presents a significant challenge for ensuring consistency between the system models and computer-aided design/engineering (CAD/CAE) models. This paper concentrates on the crucial technical issue of systematically verifying the consistency between the two models.

LPMs use a network of lumped components to represent the structure of engineering systems, such as masses, springs, and dampers, which are governed by a system of ordinary differential equations (ODEs) in terms of variables that vary with time [5,6]. On the other hand, DPMs explicitly consider the geometric and material properties of engineering systems, which are governed by a system of partial differential equations (PDEs) that take into account variables varying with both time and space [7]. Numerical discretization methods can be used to approximate these PDEs by large systems of ODEs upon approximating the spatial continuum in a finite basis [8].

Figure 1 depicts the process of designing a 3D suspension mechanism based on a system design. After creating an LPM in Modelica [9], mechanical parts are used to realize the lumped components and obtain the expected behaviors. The stiffness and damping of the absorber are modeled by a spring–damper pair, without considering the precise geometric realization. Two potential geometric realization options for the absorber are presented in the library. While there might be various options for the absorber, a designer must verify that the selected option behaves as intended by the lumped component (① in the figure) before the assembly process or design qualified parts guided by the three model consistency conditions explained in Sec. 3. Additionally, it is crucial to confirm that the final geometric assembly behaves as intended by the LPM (depicted as ② in the figure) to ensure a qualified design. Currently, the only reliable way to compare design behaviors is to simulate both the LPM and DPM and compare their differential equation solutions a posteriori, which is computationally prohibitive for large-scale models, not to mention the additional challenges related to selecting appropriate time-steps, stability, and convergence. The goal of this paper is to propose a systematic method to check consistency between the system models and CAD/CAE models.

### 1.2 Contributions.

Our main technical contributions are:

We propose a generalizable definition of

*consistency*between mechanical LPMs and DPMs that considers both model*specifications*and*solutions*, taking into account factors such as mass, initial and boundary conditions (ICs/BCs), and the behavior of interest (BoI).We develop a

*simulation-free*scheme for checking the consistency between LPMs and DPMs based on the definition proposed above. This idea is to compute a priori error bounds between the solutions of both model types by comparing their parameters, circumventing the costly process of solving differential equations.We adopt a model order reduction (MOR) technique that a priori guarantees accuracy, stability, and convergence in the simulation-free scheme to improve the efficiency of consistency analysis for large-scale models. The MOR enables the analysis of large-scale designs in a computationally efficient manner.

### 1.3 Outline.

In Sec. 2, we review the prior works on system-based geometric design, physical model solution comparison, and MOR. In Sec. 3, we define the problem of model consistency analysis. In Sec. 4, we present a simulation-free scheme to compare LPMs and DPMs. We illustrate the application of our scheme to various mechanical problems in Sec. 5. Finally, in Sec. 6, we summarize the implications of our proposed scheme and suggest potential future research directions.

## 2 Related Work

### 2.1 System-Based Geometric Design.

The system-to-geometry design strategy has been shown by Ulrich [10] to reduce complexity and clarify the problem-solving process. Most system-based geometric design approaches rely on a one-to-one correspondence between the lumped component in the LPM and the geometric part in a solid model repository. For instance, Finger and Rinderle proposed a component database that contains correspondence between ports of bond graphs and geometric parts for a specific class of mechanical designs [11]. Engelson et al. developed an integrated environment that combines geometric design and system modeling tools to assist engineers in constructing and verifying large, moving rigid-body assemblies [12]. However, these approaches may result in unrealistic designs by ignoring *function-sharing*. Ulrich’s work stands out as an exception [10], as it demonstrates a systematic approach to merge multiple functions, abstracted by different system-level components, in fewer geometric parts. Despite the numerous approaches proposed for system-based geometric design [10,12–16], none of them formally introduce the concept of consistency between system and geometric designs, nor do they provide a systematic approach for assessing the validity of the geometric design with respect to the target system behavior.

### 2.2 Solution Comparison Between the Lumped Parameter Model and the Distributed Parameter Model.

Comparing the solutions of LPMs and DPMs is a common requirement in model conversion problems, where the two models must be converted to each other with tolerable differences in their respective solutions. In computer graphics, for instance, deformable objects like cloth fabrics and soft tissues are often converted to lumped mass-spring models for faster simulations due to the simplicity and efficiency of LPMs [17–19]. Gelder [20] developed a lumped-spring element based on the geometric angle and length information of 2D linear triangular finite elements, while Baudet et al. [21] extended the method to rectangular finite elements. In both cases, the solutions of the converted LPM and the original finite element model can be directly compared. Natsupakpong and Çavuşoğlu [22] proposed a method for converting 3D deformable objects from finite element models to LPMs by minimizing the difference of the stiffness matrices of the two models. However, the dimension of the LPM solution is often much smaller than that of the DPM after spatial discretization, making direct comparison challenging. In other words, a one-to-one correspondence between the solution dimensions of the two models is typically not guaranteed. This limits the applicability of existing methods to the model solution comparison problem addressed in this paper.

### 2.3 Model Order Reduction.

The field of MOR has seen the development of various techniques [23–27] to approximate a given “full-order” model (FOM) in a numerically efficient and stable manner while preserving certain desired properties. The balanced truncation (BT) method, for instance, removes weakly controllable and observable states from the FOM [28], but can be computationally expensive for large-scale models [29]. The rational Krylov subspace (RKS) method, on the other hand, approximates the FOM transfer function by matching a few significant terms of its Taylor series expansion [30], but does not have a standard rule for choosing the frequency, which is the frequency around which the Taylor series expansion is made. To address this limitation, the iterative rational Krylov algorithm (IRKA) was developed, which iteratively updates the expansion frequency using the reflection of poles of the updated reduced order model (ROM) about the imaginary axis at each iteration step until the difference between the transfer functions of FOM and ROM is minimized [31].

However, the IRKA algorithm cannot ensure the error converges to a local minimum [32]. To overcome this limitation, the CUmulative REduction (CURE) scheme was proposed, which adaptively chooses the expansion frequency and incrementally increases the scale of the ROM by monotonically decreasing the norm of the error transfer function to zero through an accumulation process [33]. Furthermore, a stability-preserving, adaptive rational Krylov (SPARK) algorithm [33] was developed to maintain model stability and is usually embedded in the CURE scheme to generate a family of stable ROMs whose orders are increased by sequential accumulation in a single MOR process. This embedded CURE scheme with SPARK algorithm (hereafter abbreviated by SPARK + CURE) was used in Ref. [34] to generate a family of physically interpretable multi-fidelity surrogate LPMs for physical systems governed by PDEs.

The $SPARK+CURE$ method has several key advantages, namely: (1) automatic search for proper expansion frequencies, (2) preserved model stability, (3) guaranteed error convergence, and (4) a priori error bound. This method is particularly suitable for large-scale mechanical systems, as it can significantly improve the time efficiency of model consistency analysis. In Table 1, we compare the important properties of the BT, RKS, IRKA, and SPARK + CURE methods, the advantages of the latter are in bold. For our proposed simulation-free scheme (Sec. 4), we adopt the SPARK + CURE method to ensure accuracy, stability, and convergence in our model consistency analysis.

BT | RKS | IRKA | $SPARK+CURE$ | |
---|---|---|---|---|

Numerical efficiency | Low | High | High | High |

A priori error bound | Yes | No | No | Yes |

Auto order decision | No | No | No | Yes |

Maintain stability | Yes | No | Yes | Yes |

Guarantee convergence | Yes | No | No | Yes |

BT | RKS | IRKA | $SPARK+CURE$ | |
---|---|---|---|---|

Numerical efficiency | Low | High | High | High |

A priori error bound | Yes | No | No | Yes |

Auto order decision | No | No | No | Yes |

Maintain stability | Yes | No | Yes | Yes |

Guarantee convergence | Yes | No | No | Yes |

## 3 Problem Formulation

This section outlines the problem formulation for the model consistency analysis procedure, which involves two model types: the LPM, which provides a system-level description of physical systems, and the DPM, which incorporates the spatiotemporal distribution of materials. The definition of these models and the relationship between their model specifications, behaviors, and ICs/BCs are crucial for formulating the problem, and will be discussed below.

### 3.1 Definitions.

We focus on linear time-invariant (LTI) translational mechanical LPMs that can be modeled as networks of interconnected lumped components such as masses, springs, and dampers. To uniquely determine the system behavior, initial mass displacements, velocities, and source terms must be specified, which give rise to a system of ODEs in the time domain. We assume that any physically meaningful mass–spring–damper network can be realized by at least one DPM, which has a continuous material distribution and geometry embedded in 3D space. The DPM can be divided into pieces such that there is a one-to-one correspondence between the spatial integration of mass density and lumped masses. The effective stiffness and damping of the DPM correspond to those of the mass–spring–damper network, and the spatial integration of ICs/BCs and body effects of different pieces correspond to the ICs and source terms of lumped masses. The DPM behavior is described by PDEs defined over a region of spacetime. To compute the effective stiffness and damping, solving PDEs is usually unavoidable. The solution generally involves spatial discretization and numerical methods such as finite difference, element, and volume analysis [8].

Based on the above definitions, we can define *consistency* between an LPM and a DPM as the satisfaction of the following three conditions:

(C1) The total lumped mass value of the LPM matches the spatial integration of the mass density of the DPM.

(C2) The spatial integration of the ICs/BCs, as well as the body effects of the DPM, match the ICs and source terms of the LPM.

(C3) The BoI of the DPM matches the BoI of the LPM in terms of spatial integration.

In other words, for an LPM and a DPM to be consistent, the mass, ICs, source terms, and BoI of the LPM must be equivalent to the mass distribution, initial and boundary conditions, and BoI of the DPM, respectively.

### 3.2 A Mechanical Example.

Figure 2 depicts an example of an LPM and its corresponding DPM realization, where the geometry of a 3D suspension system is considered. The LPM is composed of a mass, two springs, and a damper, and is subjected to a time-varying external force $f(t)$. The network is connected to the ground via a spring. On the other hand, the DPM is an assembly of several solid parts that have linear-elastic material distributed in 3D space. The top surface of the assembly is subjected to time-varying pressure $p(x,t)$, where $x$ is the spatial 3D coordinates. The entire assembly is fixed to the ground. To ensure consistency between the LPM and DPM, we require that the following conditions are satisfied: the mass of the DPM matches the total mass of the LPM, i.e., $m=\u222b\rho (x)dV$ where $\rho (x)$ is the material density and $dV$ represents an infinitesimal volume element; the external load on the LPM matches the total force on the DPM, i.e., $f(t)=\u222bp(x,t)dS$ where $dS$ represents an infinitesimal surface element; the pre-specified displacement $u\xaf(x,t)$ of the LPM (which is fixed to the ground) matches that of the DPM, i.e., $w\xaf(t)=(1/Sc)\u222bu\xaf(x,t)dS=0$, where $Sc$ is the contact surface area between the tire and the ground; and the solved displacement of the LPM matches that of the DPM, i.e., $w(t)=(1/Sp)\u222bu(x,t)dS$, where $Sp$ is the top surface area of the DPM assembly.

### 3.3 Formulation.

To analyze consistency between LPM and DPM, one has to verify if the conditions (C1) through (C3) are satisfied.

For practical purposes, we introduce an upper-bound (i.e., inequality constraint) to relax the strict equality constraint of BoI in (C3). Noting that an LPM is a ROM of a DPM, its solution can be at best an approximation of the DPM solution, and hence, it is not reasonable to expect a strict equivalence between the solutions of LPM and DPM.

Given an LPM and a DPM with matched mass properties, ICs/BCs, and source terms, our objective is to check if the overall error between pairs of LPM BoI quantities $wi(t)$ and the corresponding DPM BoI quantities $ui(t)$ for $i=1,2,\u2026,n$ is less than an acceptable tolerance $\u03f5>0$. Symbolically, this can be represented as $\u2211i=1n\Vert ui(t)\u2212wi(t)\Vert \u2264\u03f5$. Any commonly used vector norm such as $L1$, $L2$, and $L\u221e$ can be used as the norm $\Vert \u22c5\Vert $. However, for the sake of specificity, we use $L\u221e$ to provide a measure of the maximum deviation between the DPM and LPM BoI quantities.

Below, we present a simulation-free scheme that compares the BoIs of a given pair of LPM and DPM in terms of the $L\u221e$ norm, for the purpose of model consistency analysis.

## 4 Simulation-Free Scheme

Assuming that an LPM and a DPM are provided, along with their matching mass properties, ICs/BCs, and source terms, as well as the BoI quantities’ matching form (i.e., spatial integration over the boundary surface), a scheme is presented below for calculating the a priori error between the corresponding BoI quantities of the two models, without requiring explicit numerical solution of differential equations.

We start by deriving the governing equations for both the LPM and the DPM, which are respectively described by a system of second-order ODEs and PDEs, respectively. To enable a comparison between the two models, we discretize the PDEs of the DPM using spatial discretization methods [35], thereby converting them to a system of second-order ODEs. We then transform the given ODEs of the LPM and the resulting ODEs (semi-discretized PDEs) of the DPM into their respective state-space forms through linear transformation [36], as illustrated below:

- LPM ODEs in state-space form:(1)$Elx\u02d9l(t)=Alxl(t)+Blhl(t)yl(t)=Clxl(t)$
- DPM ODEs (semi-discretized PDEs) in state-space form:(2)$Edx\u02d9d(t)=Adxd(t)+Bdhd(t)yd(t)=Cdxd(t)$

Since both systems of ODEs are LTI, we can match them using linear projections:

The source terms are related by $Bl(t)hl(t)=\Gamma nBd(t)hd(t)$, where $\Gamma n$ represents a projection matrix that maps the discrete form of forces of DPM to the lumped forces of LPM.

The initial displacement and velocity terms are related by $\Gamma Ixd(0)=xl(0)$ and $\Gamma Ix\u02d9d(0)=x\u02d9l(0)$, where $\Gamma I$ denotes the projection matrix that maps the spatially discretized ICs of the DPM to the ICs of the LPM.

The BoI terms are related by $Cd=\Gamma f$, where $\Gamma f$ is a projection matrix that maps the spatially discretized BoI of DPM to the BoI of LPM.

*a*)

*b*)

Notably, the integration term in (6) is constant for a given source term $hd(t)$, as the BC of DPM is time-invariant. This means that as the norm $\Vert Gd(s)\u2212Gl(s)\Vert H2$ approaches zero, the maximum difference between $yd(t)$ and $yl(t)$ also approaches zero.

In essence, the deviation between the BoIs of LPM and DPM in the time domain can be approximated by the deviation between $Gd(s)$ and $Gl(s)$ in the frequency domain using the $H2$ norm. An important observation is that the computation of $\Vert Gd(s)\u2212Gl(s)\Vert H2$ involves solving algebraic equations instead of ODEs, as discussed in Ref. [39]. This means that we can avoid solving ODEs altogether. To speed up the computation of these algebraic equations, we use the ROM of the DPM, as shown in the equation below.

In order to enhance the computational efficiency of solving large-scale algebraic equations arising from DPMs, a MOR method is employed. Specifically, we adopt the SPARK + CURE method from Ref. [33], which generates a small-scale surrogate model for the DPM. Due to space constraints, please refer to Panzer’s original Ph.D. thesis [33] for more details about this method. This surrogate model has the same source terms $ud(t)$ as the DPM in (2), and its transfer function is denoted by $Gr(s)$. The SPARK + CURE method guarantees an upper-bound $\Vert Gd(s)\u2212Gl(s)\Vert H2\u2264\epsilon \xaf1$ between the transfer functions of the DPM and the surrogate model. Unlike the spatial-discretized DPM, we assume that the given LPM in (1) has a small scale, so we do not apply MOR to it.

Figure 3 provides a visual representation of the five-step simulation-free scheme proposed in this study, where an LPM with masses, springs, and dampers and a DPM with 3D suspension are used as examples. The approach is outlined step-by-step as follows. Given the ODEs of the LPM and the PDEs of the DPM:

Use a spatial discretization (e.g., finite element [35]) method to convert the PDEs into a system of ODEs.

Convert the ODEs of both models to state-space form.

Use the SPARK + CURE method to reduce the scale of the DPM and ensure a priori guaranteed $H2$ error $\epsilon \xaf1$ between the DPM and its surrogate model.

Calculate the exact $H2$ error $\epsilon \xaf2$ between the surrogate DPM and the LPM.

Estimate the upper-bound of the $H2$ error $\epsilon \xaf$ between the transfer functions of the given LPM and semi-discretized DPM using the triangular inequality and verify whether it falls within the acceptable tolerance.

## 5 Applications

In this section, we demonstrate the practical application of the proposed simulation-free scheme by using it to compare the solutions of the LPM and the DPM for two mechanical designs: a bracket (Fig. 4(b)) and a frame (Fig. 4(c)). Both designs have linear isotropic material properties, and their mass properties, ICs/BCs, and the form of the BoI quantities are pre-specified to ensure consistency between the two models.

### 5.1 Model Specifications.

The same LPM topology is utilized for both mechanical designs, as shown in Fig. 4(a), but with different parameters for the lumped components. The values are provided in Table 2, using SI base units. Both lumped masses have zero initial displacements and velocities.

$m1$ | $m2$ | $k1$ | $k2$ | $r1$ | $r2$ | $f1$ | |
---|---|---|---|---|---|---|---|

LPM for bracket | $3.8465\xd7105$ | $3.512\xd7103$ | $3.316\xd7104$ | $4.688\xd7103$ | $1.4697\xd7105$ | $2.9052\xd7103$ | $0.005$ |

LPM for frame | $7.997\xd7105$ | $6.9139\xd7104$ | $4.8561\xd7108$ | $2.2308\xd7108$ | $2.8102\xd7107$ | $1.5075\xd7105$ | $2186.56$ |

$m1$ | $m2$ | $k1$ | $k2$ | $r1$ | $r2$ | $f1$ | |
---|---|---|---|---|---|---|---|

LPM for bracket | $3.8465\xd7105$ | $3.512\xd7103$ | $3.316\xd7104$ | $4.688\xd7103$ | $1.4697\xd7105$ | $2.9052\xd7103$ | $0.005$ |

LPM for frame | $7.997\xd7105$ | $6.9139\xd7104$ | $4.8561\xd7108$ | $2.2308\xd7108$ | $2.8102\xd7107$ | $1.5075\xd7105$ | $2186.56$ |

Two mechanical parts, namely, a bracket and a frame, have been designed to implement the LPM instances. The DPMs for both designs have zero initial displacements and velocities. We discretize the bracket and the frame using second-order tetrahedral finite elements (Fig. 5) to generate DPM ODEs. The top surfaces of the bracket and the frame are subjected to pressures that vary with time. The time-varying characteristics of the pressures are depicted in Figs. 6(a) and 6(b) for the bracket and the frame, respectively.

In order to compare the LPM and DPM models of the two designs, it is necessary to ensure that their mass properties, ICs/BCs, and BoI quantities match. The requirements for this matching are explained in (C1)–(C3) in Sec. 4, and an example of the matching process is provided in Sec. 3.2. To maintain consistency, we will follow the same matching procedure as the example in Sec. 3.2.

### 5.2 Surrogate Distributed Parameter Model of the Bracket.

Here, we demonstrate the effectiveness of our simulation-free approach in enhancing the efficiency of model comparison for large-scale models. We used tetrahedral finite elements to discretize the bracket’s geometry at three different resolutions (only one sketch is shown in Fig. 7(a)), resulting in three sets of state-space equations with a state variable count ranging from approximately 3000 to 6000. We then applied the SPARK + CURE method to generate three sets of ROMs that have a priori relative $H2$ error bounds just below 1%.

To assess the quality of the ROMs, we used the backward Euler method [7] to simulate the DPM and all its surrogate models, using the same solver and time-steps. The simulation results are presented in Figs. 7(a) through 7(c), where the abscissa denotes computation time and the ordinate denotes the average displacement of the bracket’s top surface. Due to limited space, only a few ROM curves are labeled. The solver employed in the simulations was the “sparse state space” toolbox in Matlab [40], which can analyze high-dimensional dynamical systems with state-space dimensions of $O(104)$ or more. This solver preserves system matrix sparsity, allowing it to use computationally efficient operations, such as sparse lower–upper decompositions [41], that would otherwise be infeasible or time-consuming. It is important to note that the simulations were only used to illustrate the accuracy of the surrogate models and were not necessary for our proposed model comparison method.

### 5.3 Error Analysis.

As the order of ROM increased, which corresponds to an increase in the number of state variables of ROM, the simulated response of the ROM more closely approached that of the DPM for all cases. To quantify the error introduced by the MOR process, we calculated two different measures, as shown in Fig. 8: the relative $H2$ error of the transfer functions and the root mean square error (RMSE) of the simulated response, whose formulas can be found in Ref. [42]. The SPARK + CURE method guaranteed the strictly monotonic decay of the bound of the a priori relative $H2$ error, as shown in Figs. 8(d)–8(f). However, it did not provide a rigorous a priori bound on RMSE, as shown in Figs. 8(a)–8(c), although an overall reduction pattern was commonly observed.

Given that the pre-conditions ensure the matches of masses, ICs, and BCs, the LPM and DPM would be considered consistent if the upper-bound of error $\epsilon \xaf=0.0463$ is within the user-defined tolerance. To compare the numerical solutions between the LPM and DPM, we plot the errors in Fig. 9(b), while the numerical solutions are compared in Fig. 9(a). The maximum difference between the model solutions over the entire time domain is less than $4.8197\xd710\u22128$, which is one order of magnitude smaller than the maximum model output $O(10\u22126$).

### 5.4 Time Efficiency Analysis.

Through simulations of the DPM and all ROMs for three cases of the bracket, we have obtained two polynomial curves (Fig. 7(d)) that fit the relation between computation time and the order of DPM for both a posterior solution comparison based on direct ODE simulations and our simulation-free scheme. The curves intersect at a certain point, indicating that if the order of the discretized DPM is lower than 4637, comparing solutions through simulations is more time-efficient. However, if the order is higher than 4637, our simulation-free scheme will be more time-efficient, with an accuracy loss of at most 1%. It should be noted that if a smaller user-selected relative $H2$ error bound is used (i.e., $<1%$), a new cross-over point that slightly shifts toward a larger order of DPM will be found.

### 5.5 Surrogate Distributed Parameter Model of the Frame and Error Analysis.

The agreement between the LPM and DPM can be determined by the user-defined threshold for the relative error. When the threshold is greater than $0.8229$, the models are in agreement. If the threshold is smaller than $0.8229$, the models are not in agreement. In Fig. 11(a), we compare the simulation results of these two models, and their numerical solutions difference is depicted in Fig. 11(b). It can be observed that the difference is significantly smaller than the displacement of design models in most of the time-steps.

## 6 Conclusion

The paper introduces a model consistency concept and presents a simulation-free scheme for comparing the behavior of the LPM and the DPM. By leveraging the SPARK + CURE MOR method, the proposed approach can address the computational challenges associated with large-scale models. Our analysis of two mechanical models demonstrates that our simulation-free scheme provides relatively tight $H2$ a priori error bounds for behavior comparison. Furthermore, we show that our proposed scheme is more time-efficient than direct numerical simulation if the order of the discretized DPM is above a certain threshold, with increasing efficiency as the DPM order increases. The proposed approach can be applied to various engineering applications to enhance the system-based geometric design and facilitate integration between system-level and geometric designs.

However, we acknowledge that our framework has some limitations. One limitation is that the scheme to compare model solutions does not account for errors caused by spatial discretization of DPMs, and users must select an appropriate spatial discretization method and resolution. Additionally, the CURE scheme can only be applied to reduce the order of dissipative models, and the a priori guarantees provided by SPARK + CURE are only valid for LTI models. There are two ways to extend our proposed scheme to nonlinear models. The first approach is to divide the nonlinear MOR models into piecewise linear MOR models, and then apply the SPARK + CURE method, similar to the trajectory piecewise linear method [43]. The second approach is to use MOR methods designed specifically for nonlinear models, such as the symplectic MOR methods [44], the MOR methods based on the proper orthogonal decomposition [45], the Gaussian process regression [46], etc. However, many of these methods do not offer rigorous a priori error guarantees.

## Acknowledgment

This material is based upon work supported by the Defense Advanced Research Projects Agency (DARPA) under Agreement Nos. HR00111990029 and HR00112090065. The responsibility for errors and omissions lies solely with the authors.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.