Hybridsimulationisamethodusedtoinvestigatethedynamicresponseofasystem subjectedtoarealisticloadingscenariobycombiningnumericalandphysicalsubstruc- tures. To ensure high ﬁdelity of the simulation results, it is often necessary to conduct hybrid simulation in real-time. One of the challenges arising in real-time hybrid sim- ulation originates from high-dimensional nonlinear numerical substructures and, in particular, from the computational cost linked to the computation of their dynamic re- sponses with suﬃcient accuracy. Itis oftenthe case thatthe simulation time-step must bedecreasedtocapturethedynamicbehaviorofnumericalsubstructures,thusresulting inlongercomputation. Whensuchcomputationtakeslongerthantheactualsimulation time, time delays are introduced and the simulation timescale becomes distorted. In suchacase,theonlyviablesolutionfordoinghybridsimulationinreal-timeistoreduce the order of suchcomplex numericalsubstructures. In this study, a model order reduction framework is proposed for real-time hybrid simulation, based on polynomial chaos expansion and feedforward neural networks. A parametric case study encompassing a virtual hybrid model is used to validate the framework. Selected numerical substructures are substituted with their respective reduced-order models. To determine the robustness of the framework, parameter sets are deﬁned to cover the design space of interest. A comparison between the full- and reduced-orderhybridmodelresponseisdelivered. Theattainedresultsdemonstratethe performance of the proposed framework. Keywords: real-time hybrid simulation, model order reduction, polynomial chaos expansion, feedforwardneural networks,dynamicresponse. 1. Introduction Hybrid simulation (HS) is a method used to investigate the dynamic response of a system subjected to a realistic loading scenario [1]. Across diﬀerent engineering DRAFT disciplines,HSisalsoknownashardware-in-the-loop(HiL)[2]orsystem-in-the-loop. It is one of the methods used for Model-based Development [3]. The system-under- consideration consists of multiple individual sub-components, out of which one or more are selected to be tested physically, and are called the physical substructures (PS), while the remaining components, called the numerical substructures (NS), are simulated numerically. The coupled PS and NS constitute the so-called hybrid model. Inmoredetail,thehybridmodeldynamicresponseisevaluatedonline,basedonastep- by-step numerical solution of the governing equations of motion of the hybrid model, combined with experimental responses obtained from the PS. The coupling of NS and PS is achieved by an arrangement of transfer systems, such as servo-controlled motors and actuators. These transfer systems are responsible for synchronizing the dynamic boundary conditions of the substructures at their interfaces throughout the entire HS. TheadvantageofHSliesinthefactthatitcanbeusedtoexaminetheinnerworkingsof asystem’ssub-component,beyonditslinearregime,withouttestingtheentirefull-scale system. Asaconsequence,realisticresultsfromthetestedcomponentareobtainedata much lower experimentation cost. The report by Schellenberg and Mahin [4] provides a comprehensivereview ofHS methodsand algorithms. Frequently, HS is performed on a distorted timescale, e.g. by slowing down the rate of PS testing to accommodate the power or speed of the transfer system in use. These HSs correspond to the so-called pseudodynamic tests [5]. However, in order to increasetheﬁdelityandtrustinHSoutcomes,HSshouldbeconductedasclosetoreal- time as possible. This is especially important for loading-rate-sensitive substructures. Even though real-time hybrid simulation (RTHS) has many advantages, it comes with numerous challenges. The ﬁrst challenge originates from the inherent dynamics of the transfersystemexcitingthePS,asitintroducestimedelaysintheoverallHSthatalterthe dynamic response of the tested hybrid model. This problem is addressed by designing control systems to adequately compensate for these delays. Several control strategies for RTHS have been proposed and are available in the literature. Some of the recently developed novel approaches can be found in [6, 7,8,9]. Nonetheless, challenges in RTHS do not arise only from the transfer system of the PS but also from the NS side, such as from the computational power needed to compute the NS responses. In order to capture the dynamic behavior of interest of a high-dimensional nonlinear numerical model, it is often necessary to use small time-steps for the numerical solver. However, the smaller the time-step of the simulation, the larger the computational power needed to compute it. In RTHS, when the required computation time becomes longer than the actualsimulationtime,timedelaysareintroducedtothehybridmodelriskingnotonly todecreasetheﬁdelityoftheHSresultsbutalsotodrivetheoverallHSintoinstability. Therefore, insuch cases,the only viablesolution is toreduce theorder of theNS tobe able toconduct theHS inreal-time. Several model order reduction (MOR) techniques can be found in the literature. Properorthogonaldecomposition(POD)[10, 11,12],also knownas principalcompo- nentanalysis(PCA),anditsdiscreteversion,namelythesingularvaluedecomposition (SVD), are mostly combined with the method of snapshots [13, 14] and have been extensivelyusedforextractionofmodeshapesfromhigh-dimensionalsystems,among other techniques. A technique similar to POD, also using the method of snapshots but following a data-based instead of a model-based approach, is the dynamic mode 2 DRAFT decomposition(DMD)[15]. Inparticular,themaindiﬀerencebetweenPODandDMD isthattheformerreliesuponatime-averagedspatialcorrelationmatrixwhilethelatter representsthetemporaldynamicswithhigh-degreepolynomials[16]. Anadvancement of POD was also proposed based on the use of separated representations, the so-called proper generalized decomposition (PGD) [17]. More recently, the component-mode synthesis (CMS) approach [18, 19] has been used for non classically damped linear systems [20] as well as in HS as a method to reduce the order of the comprising NS andPS[21]. Furthermore,aquadraticmanifold[22]andanon-intrusive[23]approach were proposed for MOR of nonlinear structural systems. MOR techniques have also been used for Bayesian ﬁnite element model updating [24] and for reliability-based design problems[25]. In this study, a MOR framework originating from a data-driven regression is pro- posed to reduce the order of high-dimensional nonlinear NS in HS. Two diﬀerent methodologies for MOR are addressed; a polynomial chaos expansion (PCE) and a feedforward neural network (FFNN). In both cases we treat the NS, whose order we wanttoreduce,asa blackbox,meaningthatnopriorknowledgeofitsfull-ordermodel (FOM) dynamics is required. Instead, we map the outputs of the respective NS to its inputsbysimpler-to-evaluatefunctionsthanthoseoftheoriginalFOM.Thegoalofthe MORframeworkistoadequatelycapturethedynamicbehaviorofhigh-dimensionalNS at a much lower computational cost and in shorter computation time, and thus enable real-time HS. A parametric case study encompassing a virtual hybrid model is em- ployed to validate the proposed MOR framework. Selected NS are replaced with their respective FFNN- and PCE-based reduced-order models, while speciﬁc parameters of the hybrid model are varied to investigate the framework’s robustness to diﬀerent hy- brid model conﬁgurations. The corresponding reduced-order hybrid model responses are compared to the ones from the FOM. Results demonstrate the eﬀectiveness of the proposed MORframework. The paper is organized as follows. Section 2introduces the MOR framework describingthePCEandFFNN.Section 3presentsthecasestudywiththevirtualhybrid model used to validate the framework. Section 4discusses the obtained results and Section5presentsthe overall conclusionsof thisstudy. 2. Polynomial chaos expansion and neural networks as model order reduction techniques Considering an input vector Xand a computational model Y=M(X ), data- driven regression algorithms formulate a map M:X/mapsto→Ybased on an obtained samplesetofinputpoints X=/braceleftbig x(1),...,x(N)/bracerightbigTandoftherespectiveoutputvalues, i.e., model evaluations, Y=/braceleftbig y(1),...,y(N)/bracerightbigT. The set ofX,Yrealizations corre- sponds to the so-called experimental design (ED). Regression techniques encompass linearregression,kernelmethods,neuralnetworksandgraphicalmodelsamongothers. An overview of these techniques can be found in [26, 27]. As mentioned above, the proposed MOR framework of this study explores use of PCE and FFNN to this end, with thesemethods elaboratedin whatfollows. 3 DRAFT 2.1. Polynomial chaos expansion PCE is a well-known uncertainty quantiﬁcation spectral method used to substitute the dynamics of an expensive-to-compute numerical model with an inexpensive-to- compute surrogate (a.k.a. metamodel), representing the outputs of the model by a polynomial function of its inputs [28, 29]. It is proven to be a powerful surrogate technique used in a wide variety of engineering contexts to replicate the dynamic response of complex high-dimensional models [30, 31,32,33]. Hence, it can be seen asapromisingtechniqueforMORofNSinHS.Furthermore,tothebestoftheauthors’ knowledge,suchuse ofPCE withinthe HScontext isnew. In more detail, given a random input vector X∈RMwith statistically inde- pendent components expressed by the joint probability density function (PDF) fX and a ﬁnite variance computational model Y=M(X ), withY∈Rsuch that E[Y2] =/integraltext DXM2(x)fX(x)dx<∞,the PCEofM(X )follows: Y=M(X )≈MPCE(X) =/summationdisplay α∈NMyαΨα(X). (1) The PCE function is built on the Ψα(X)multivariate orthonormal polynomial basis with respect to the input vector fX. The degree of the Ψαpolynomials components is identiﬁed by the α= (α 1,...,αM),α∈NMmulti-index for each of the input variables,while yαcorrespondsto thepolynomialcoeﬃcients. 2.1.1. Polynomial basis Themultivariatepolynomialsareconstructedasatensorproductoftheirunivariate orthonormalpolynomials φ(i) k(xi)suchthat Ψα(X) =/producttextM i=1φ(i) αi(xi). Thelattermeet the/angbracketleftBig φ(i) j(xi),φ(i) k(xi)/angbracketrightBig =/integraltext DXiφ(i) j(xi)φ(i) k(xi)fXi(xi)dxi=δjkorthonormality criteria, where icorresponds to the input variable, jandkto the polynomial degree, fXi(xi)totheith-inputmarginalPDFand δjktotheKroneckersymbol. Theselection of the univariate orthonormal polynomial families depends on the marginal PDF of each input variable to which they are orthogonal, e.g., if an input variable follows the uniform/ Gaussian distribution then the Legendre/ Hermite orthogonal polynomial familyis usedrespectivelyfor thisspeciﬁc input variable[29, 34]. 2.1.2. Truncation schemes Once the univariate polynomials families are selected for each input variable, the nextstepistheconstructionofthePCEfollowingEq.(1). However,becausetheterms of the sum in Eq. (1) are inﬁnite, it is often truncated to a ﬁnite number of terms for practical reasons. Hence the truncated basis is deﬁned as A⊂NMand the PCE of M(X )admits: MPCE(X) =/summationdisplay α∈AyαΨα(X). (2) The performance of PCE is closely connected with the chosen truncation scheme. Underﬁttingoroverﬁttingispossiblewhentoomanytermsarediscardedorintroduced, respectively [32]. In the standard basis truncation scheme [29] the maximum degree 4 DRAFT p∈N+is deﬁned, such that the degree of each polynomial is capped by this value. Therefore,thestandardbasistruncationschemeconsistof/parenleftbigM+p p/parenrightbig elementsandfollows: AM,p=/braceleftbig α∈NM:|α|≤p/bracerightbig , (3) whereMcorresponds to the number of input variables of M(X )and|α|=/summationtextM i=1αi to the total degree of all polynomials in Ψα. Alternate truncation schemes have been developed namely, the maximum interaction and hyperbolic truncation [35], and are usedbasedonthescopeofeachapplication. Inthisstudythestandardbasistruncation schemeof Eq.(3) isemployed. 2.1.3. PCE coeﬃcientcalculation Once the truncation scheme is determined, the next step is the calculation of coef- ﬁcientsyαin Eq. (2). Various methods have been developed to do this [34]. In the proposed framework, the least-squares minimization technique [36] is used. This is an non-intrusivetechnique,meaningthatthecoeﬃcientsareobtainedafterpost-processing the EDpoints. More speciﬁcally, Eq.(2) isreformedas: MPCE(X) =P−1/summationdisplay j=0yjΨj(X) +εP=yTΨ(X ) +εP, (4) whereΨ(X ) ={Ψ0(X),..., ΨP−1(X)}Tare the multivariate orthonormal poly- nomials,y={y0,...,yP−1}Tthe PCE coeﬃcients, εPthe truncation error, and P=/parenleftbigM+p p/parenrightbig . Then, thePCE coeﬃcientsare obtainedby: ˆy= arg min E/bracketleftBig/parenleftbig yTΨ(X )−M(X )/parenrightbig2/bracketrightBig . (5) The solutionof Eq.(5) isattained byOrdinaryLeast-Squares(OLS) andfollows: ˆy=/parenleftbig ATA/parenrightbig−1ATY, (6) where the experimental matrix is Aij= Ψj(x(i))withi= 1,...,N andj= 0,...,P−1. 2.2. Feedforward neural network Neural networks are an extremely versatile group of methods used both for regres- sion and classiﬁcation problems. The most straightforward and widely used variety of neural networks are the fully connected FFNN. Such a FFNN consists of layers of "neurons",wherethevalueoftheneuronsineachlayeriscalculatedasamatrixmulti- plicationbetweenaweightmatrixandthevectorofneuronvaluesatthepreviouslayer, followedbytheapplicationofanelement-wisenonlinear"activation"function. FFNNs have proven highly capable at approximating a broad range of functions and can be showntobe,inthelimitofinﬁniteneurons,universalapproximators[37, 38]. Assuch, they are an obvious candidate for metamodelling, and have been widely demonstrated in many such contexts [39, 40,41]. Much of the popularity of FFNNs is due to their scalability to very large and high dimensional datasets, and, especially in the scope 5 DRAFT of deep learning methods, power for dealing with highly nonlinear relations [42, 43]. This also makes the FFNNs a very attractive option for metamodelling in a hybrid testing context, in which they have been demonstrated in limited number of cases to date [44,45]. 2.2.1. NetworkArchitecture Much of the power of neural networks comes from the ease with which their architecturescanbeadaptedtodealwithsystemsofdiﬀerentsizesandcomplexity[46]. The basic architecture of a 3 layer FFNN is demonstrated in Figure 1. Such a network consists of the input layer, where the neurons take the values of the input vector X. Thenumberofneuronsinthislayerisnaturallyequaltothedimensionalityoftheinput vector. The ﬁnal layer of the network is the output layer, at which the neuron values are the elements of the output vector Y. The number of neurons in the output layer is equal to the dimensionality of the output of the network. Between these two layers lies a hidden layer which can have an arbitrary number of neurons, the selection of which greatly aﬀects the performance of the network. The hidden layer neurons are characterized by the activation values, which are represented in the vector A. It is possible and often desirable to use multiple hidden layers between the input and the outputlayers,creatingaso-calleddeepneuralnetwork. Deepneuralnetworksareeven more powerful for representing highly complex functions [47], however, only a single hidden layeris usedin thiswork sincesuch aFFNN isdeemed adequate. The single-hidden-layer FFNN diagram in Figure 1is mathematically described in Eqs.(7)and(8). Eq.(7)representsthecalculationoftheactivationvaluesinthehidden layeroftheFFNN,theinputvector Xismatrixmultipliedbytheweightmatrix W1and the bias vector b1is added. The activation vector Ais then calculated by the element- wiseapplicationoftheactivationfunction g1. InEq.(8),theoutputvectoroftheFFNN is calculated in a similar fashion, via matrix multiplication of the activation vector A with a second weight matrix W2, to which a second bias vector b2is added before the element-wise application of a second activation function g2. In these equations, bothweightmatrices W1andW2,andbothbiasvectors, b1andb2,arepopulatedby parameterswhicharelearnedduringthetrainingoftheFFNN.Theactivationfunctions g1andg2, on the other hand, are hyperparameters and must be chosen by the network designer. Common candidates for these activation functions include the tanh,RELU andlinearactivation functions,each withdiﬀerent associatedproperties[48]. A=g1(W1X+b1) (7) Y=g2(W2A+b2) (8) 2.2.2. NetworkTraining When using a FFNN to create a metamodel, the ﬁrst step involves generating a dataset with which to train, validate and test the network. Such a dataset consists of input and output vector pairs (points) for the system being metamodelled. For traininganeuralnetwork,thisdatasetisthenusuallyfurtherbrokendownintotraining, validation and testing datasets. The training dataset contains those input-output pairs used for optimizing the learned parameters of the system in the weight matrices and 6 DRAFT Input Layer Hidden Layer Output LayerXA YFigure 1. Architecture of a 3layer FFNN bias vectors. The validation dataset is a set of held-back input-output pairs which are used to compare FFNNs created with diﬀering hyperparameters, such as hidden layer sizeoractivationfunctions. ThisyieldsFFNNswithbetterperformanceonthetraining dataset. The use of a separate validation dataset informs on the ability of the trained FFNNtogeneralize,i.g.,demonstratewhethertheFFNNperformanceismaintainedfor input-outputpairsoutsideofthetrainingdataset,aswellasindicatewhethertheFFNN is "overﬁtting" to the training points [46]. The testing dataset ﬁnally, comprises those pointsonwhichtheﬁnalperformanceofthetrainednetworkofthechosenarchitecture is tested. In this work, the neural network was created and trained using the Matlab deep learning toolbox [49]. Since, in the context of metamodelling, the problem is one of regression, i.e., predicting a continuous output variable, the loss function of the networktobeminimizedwasselectedasthemeansquarederrorvalue,asdemonstrated in Eq. (9). In this equation, the loss value Jis calculated as the mean value of the squarederrorbetweenthetrueoutputvalue Yiforthetrainingpoint iandthepredicted output value from the network ˆYifor allNtraining data points. The network weights were trained using the Bayesian regularization algorithm. The latter is more resistant tooverﬁttingandtendstotrainnetworkswithgreatergeneralization[50]. Avalidation dataset was used to select the ﬁnal architecture of the network as that which gave the best validation performance. The ﬁnal network architecture selected for this problem wasa3layernetworkwithasinglehiddenlayer. Theactivationfunctionsforthehidden and output layersare tanhandlinearfunctions, respectively. J=1 NN/summationdisplay i=1(Yi−ˆYi)2(9) 3. Case study Given that the focus of this work is a MOR framework for RTHS, the case study employs a virtual hybrid model to conduct virtual RTHSs (vRTHSs). All of the components of the virtual hybrid model are implemented numerically. Thus, virtual PSs (vPSs) are used in place of physical specimens. In this regard, without loss of 7 DRAFT Figure 2. Hybrid modelblock diagram. generality and for the sake of clarity, transfer systems that would be required between NS and PS in a physical HSare omittedfrom thiscase study. Thevirtualhybridmodelusedinthiscasestudyrepresentsaprototypemotorcycle. This model consists of one vPS, the electronically-controlled continuously variable transmission(eCVT)ofthemotorcycle,andfourNS:i)theenginedynamicsmodel,ii) themotorcyclebodydynamicsmodelincludingthetiresandtheroad,iii)therearwheel brakingsystems,andiv)thefrontwheelbrakingsystem. Theinterconnectionsbetween NS and the vPS are illustrated in Figure 2. All substructures of the hybrid model are developed in the Simcenter Amesim software and are brieﬂy introduced below. For a more detailed description of the aforementioned models, as well as the underlying equationsgoverningeach substructure,the readeris encouraged toconsult [51]. In particular, the eCVT vPS consist of a multi-input-multi-output (MIMO) model with two sets of inputs/outputs. The ﬁrst set is connected to the engine dynamics model NS and the second to the motorcycle body dynamics NS. The latter connection corresponds to the transmission output shaft of the motorcycle. The engine dynamics model NS aims to simulate the behavior of the motorcycle combustion engine. It is a multi-input-single-output (MISO) model, the inputs being the angular velocity of the engine crankshaft ωenand thethrottle level thr, andthe output beingthe torque ofthe engine crankshaft τen. The parameters of this model are the engine moment of inertia and coeﬃcient of viscous friction. The motorcycle body NS addresses the inner body dynamics of the motorcycle along with the dynamics of the suspension and the tires as well as the road proﬁle and the environmental driving conditions. The parameters of this model include the motorcycle mass, as well as damping and stiﬀness of the front and rear tires and suspensions. The motorcycle body NS is a MIMO model with 3 sets of inputs/outputs. The ﬁrst set is connected to the eCVT vPS with the input and the output being the torque τvdand the angular velocity ωvdof the transmission outputshaft,respectively. Thesecondandthirdsetsareconnectedtotherearandfront wheelbrakingsystemNS,respectively. BothbrakingsystemsarerepresentedasMISO models. TheinputsoftherearwheelbrakingsystemNSaretheangularvelocityofthe rear wheelωrwand the applied force in the brake pedal Fbrrw, while the output is the 8 DRAFT generatedbrakingtorqueintherearwheel τrw. Similarly,theinputsforthefrontwheel brakingsystemNSaretheangularvelocityofthefrontwheel ωfwandtheappliedforce onthefrontbrakelever Fbrfw,whiletheoutputisthebrakingtorqueinthefrontwheel τfw. 3.1. Motorcycle testing scenario The performance of the motorcycle prototype is examined by testing its virtual hybrid model under predeﬁned driving, road and wind scenarios. These tests are, in fact,virtualreal-timehybridsimulations. ASimcenterreal-timesimulationplatformis used for data exchangebetweensubstructures andfor runningthevRTHSs. The case study involves a 45 sec long motorcycle driving scenario on a road with a deﬁned proﬁle and in deﬁned wind conditions. The driving scenario is described by the variation of the throttle thrand the braking forces Fbrrw=Fbrfw(front and rearwheelbrakingforcesareassumedtobeequal)describedbyEq.(10)andEq.(11), respectively, and shown in Figure 3. Note that the braking forces are measured in Newtonandthemaximumappliedthrottleis0.5,correspondingtoahalf-openthrottle. The road proﬁle, i.e. height variation of the road, is deﬁned by the sinusoidal function h(x)inEq.(12),where xdenotesthecurrentpositionofthemotorcycleinmeters. The ambientwind velocityis assumedto bezero. Each case study motorcycle driving test involves "driving" the motorcycle virtual hybridmodelfollowingthepredeﬁneddrivingscenarioalongthegivenroadinthegiven wind conditions. The position of the motorcycle is computed by solving its equation of motion. In this case study, the RK4 (fourth-order Runge–Kutta) method with a ﬁxed time-step of 0.1 msec was used to numerically solve this equation of motion. In addition to the position of the motorcycle, a number of indicative performance parameters are computed and monitored. These include global response parameters, such as the motorcycle velocity v, as well as the input and output parameters of the NSs,i.e. theangular velocitiesandtorques shownin Figure 2. thr(t) =  0.125t ,5≤t<9 0.5 ,9≤t≤13 −0.125t ,13<t≤17 0 , elsewhere(10) Fbrrw(t) =  10t,20≤t<25 50,25≤t≤32 −10t,32<t≤37 0, elsewhere(11) h(x) =/braceleftBigg 0 ,0≤x≤2 0.02/parenleftBig cos/parenleftbig2πx 3/parenrightbig −1/parenrightBig , elsewhere(12) 9 DRAFT 0 10 20 30 4000.10.20.30.40.5 01020304050Figure 3. Case study drivingscenario. 3.2. Reduced-ordernumerical substructures To demonstrate the proposed MOR framework, two of the four NS in the virtual hybrid model of the motorcycle, speciﬁcally the rear and front wheel braking systems, werereplacedwiththeirrespectivereduced-ordermodels,applyingthePCEandFFNN techniques addressed in Section 2. The braking systems were selected because they include the brake hydraulic systems and friction pads. Adequately emulating the dynamic behavior of such systems using numerical models requires small integration time-steps, making them ideal candidates for MOR. The reduced-order models of the twobraking system NSswill bereferredto asreduced-order NS(RO-NS) hereafter. To train the front and rear brake system PCE and FFNN RO-NSs, two sets of 200 points were generated for the input variables listed in Table 1, assuming the inputs are statistically independent, using the Latin hypercube sampling (LHS) method [52]. For each of the 200 generated points, full-order models (FOMs) of the front and rear braking systems were used to (accurately but expensively) compute the resulting front and rear brake torques τrwandτfw, respectively. The obtained experimental designs (EDs)for thefront andthe rearbraking systemsare: Xfw=/braceleftbig Fbrfw,ωfw/bracerightbigT Yfw={τfw}(13) Xrw={Fbrrw,ωrw}T Yrw={τrw}(14) Two types of RO-NSs were developed for each braking system, one featuring a single PCE and the other being a FFNN. The same ED was used to train, validate and testthePCEandtheFFNNRO-NSsforeachbrakingsystem,asdiscussedinSection 2. 3.3. Virtualhybrid models Fourdiﬀerentvirtualhybridmodelsofthemotorcyclewereconstructedinthecase study, diﬀering in the selection of braking system NSs. The base-line virtual hybrid 10 DRAFT Input Distribution Mean value Standard deviation Units FbrfwUniform 25 14.43 N ωfw Uniform 245 141.45 rpm FbrrwUniform 25 14.43 N ωrw Uniform 235 135.68 rpm Table 1. Marginal PDFs and their characteristics per input variable of the front and rear brake system RO-NSs. Used to constructthe EDs of Eq. (14) and Eq. (13). model features the FOM NSs of both the front and the rear braking systems. The remainingthreevirtualhybridmodelshaveRO-NSs. InthePCEvirtualhybridmodel, both braking systems are represented using PCE RO-NSs. Analogously, the FFNN virtual hybrid model uses FFNN RO-NSs. Finally, the mixed virtual hybrid model of themotorcycleusestheFFNNRO-NSforthefrontbrakesandthePCERO-NSsforthe rear brakes. Theperformanceofthethusgeneratedmotorcyclevirtualhybridmodelswastested statistically. To this end, the 11 parameters of the motorcycle model, listed in Table 2, were assigned uniform distributions with parameters identiﬁed from [53, 54,55] to represent realistic variations of these mechanical motorcycle properties. A total of 30 diﬀerent samples of the 11 motorcycle parameters were generated using the LHS method, essentially resulting in 30 diﬀerent motorcycle prototypes. These prototypes were modeled using four virtual hybrid models described above. Using the driving scenariodescribedabove,atotalof120motorcyclemodelresponsetimehistorieswere computed inthis casestudy. The variability of motorcycle prototype response in the selected driving scenario, computedusingthe(accuratebutexpensive)FOMvirtualhybridmodel,isdemonstrated in Figures 4. The 30 time histories of two indicative virtual hybrid model response variables, the front wheel brake torque τfwand the motorcycle velocity v= ˙x, are compared to their mean value. Evidently, the 30 generated motorcycle prototypes adequately represent the underlying probability space deﬁned by the parameters of Table2. All results presented hereafter correspond to the mean values of 30 virtual hybridmodel responsesimulations. Additionally, the time history responses of Figures 4can be intuitively interpreted. Figure4a depicts the time history of the front wheel brake torque τfw, generated by the employed driving scenario. As can be appreciated from Figure 3and Eq. (11), themotorcyclebrakesareactivatedbetweenthe20thand37thsecondofthesimulation and therefore the braking torque is non-zero only inside this time interval. Figure 4b illustratestheevolutionofthemotorcycle’svelocityresponse vusingtheaforementioned driving scenario. According to Figure 3and Eq. (10), for the ﬁrst ﬁve seconds of the simulation, no throttle is applied to the motorcycle and hence its velocity is zero. Betweenthe5thand17thsecondofthesimulation,thedrivingscenarioinvolvesstepping onthethrottleandthereforethemotorcyclestartstoaccelerate. Asaresult,itsvelocity isincreasing. FollowingEq.(11),asmentionedbefore,afterthe20thsecond,thedriving scenario involves pressing on the brakes and thus the motorcycle’s velocity begins to decrease andit isbrought toalmosta fullstopin the45 secondsof thesimulation. 11 DRAFT Parameter Description DistributionMean valueStandard deviationUnits KrtVertical stiﬀness rear tire Uniform 58570 11714 N/m ZrtVertical damping rear tire Uniform 11650 3495 Ns/m KftVertical stiﬀness front tire Uniform 25000 5000 N/m ZftVertical damping front tire Uniform 2134 640.2 Ns/m KrsStiﬀness rear suspension Uniform 125000 25000 N/m ZrsDamping rear suspension Uniform 10000 3000 Ns/m KfsStiﬀness front suspension Uniform 19000 3800 N/m ZfsDamping front suspension Uniform 1250 375 Ns/m MMotorcycle massUniform 300 6 Kg JEngine moment of inertia Uniform 0.0115 0.0023 Kgm2 µEngine coeﬃcient of viscous friction Uniform 0.001 0.00005 Nm/(rev/min) Table 2.Varying parametersof thehybrid model and their characteristics. 0 10 20 30 40051015202530 (a)Front wheel braking torque. 0 10 20 30 4005101520253035 (b)Motorcycle velocity. Figure 4. ThirtyFOMvirtualhybridmodelresponsestimehistoriesandtheirrespectivemeanvaluesforthe front wheel brakingtorque (a)and the motorcycle velocity (b). 12 DRAFT 4. Results and discussion The performance of the motorcycle prototype wasevaluated in tests involving a 45 sec long driving scenario. Real-time hybrid simulations lasted approximately 42 sec when the FOM virtual hybrid model was used, while they took approximately 19 sec, 20secand21secfortheFFNN,PCEandmixed(FFNN&PCE)virtualhybridmodels, respectively. Thus, by reducing the order of two braking system NSs, the computation time was halved. This is notable because the computation time of the FOM virtual hybridmodelisclosetothedurationofthedrivingscenario: itisconceivablethatifthe drivingscenariowasmorecomplex,theFOMvirtualhybridmodelwouldnotcomplete its calculations in real-time, risking time delays and distortion of the HS timescale. Therefore, reducing the order of the critical NSs using the techniques discussed in Section 2 is justiﬁed. TofurtherassesstheperformanceoftheMORframework,thenormalizedrootmean squareerror(NRMSE)isusedtocomparetheselectedhybridmodeldynamicresponse parameters. Note that the compared responses were computed from the four sets of 30 driving scenario simulations using four virtual hybrid models of the motorcycle prototype described above. First,theRO-NSsofthefrontandrearbrakingsystemsarecomparedtothefull-order model(FOM)prototype. RecallthatthesameEDs(Eqs.(13)and(14)),comprising200 input-output data points generated by the FOM of the front and rear braking systems, were used to train the RO-NSs. To assess the RO-NSs performance, the quantities that characterize the braking systems, namely the rear τrwand front wheel braking torque τfw, are compared to the same outputs of the FOM NS. Table 3displays the NRMSE between FOM NS and RO-NSs for τrwandτfw. The NRMSE is small, less than 3% fortherealwheeltorqueandlessthan7%forthefrontwheel,withoutaclearadvantage of the FFNN or the PCE RO-NSs. The NRMSE error diﬀerence likely originates with the motorcycle dynamics. Figures 5a-b present the mean time history response of the braking torques obtained from the FOM, FFNN and PCE virtual hybrid models. The comparison indicates that the responses from the reduced-order virtual hybrid models are similar and only slightly diﬀerent from the response of the FOM virtual hybrid model. Second, the response of the four virtual hybrid models is compared by examining theNRMSEsandthetimehistoriesoftheindicativemodelresponsequantities,namely the angular velocities and torques that couple the substructures of the hybrid models (ωrw,ωfw,ωen,ωvd,τenandτvd) and the motorcycle velocity v. The variables that couple the substructures of the hybrid models were investigated since deviations from their nominal values, would risk the HS outcomes and could drive the HS loop into instability. In addition, for non-virtual RTHS that include real test benches and PS, large deviations could damage the corresponding physical specimen and laboratory equipment. On the contrary, the motorcycle velocity vis chosen to be examined as it corresponds to a global response quantity that informs about the overall dynamic behavior of the motorcycle during the driving scenario. Table 4provides an overview of the NRMSE and Figures 6,7show the response quantity time histories for the four virtual hybrid models. Evidently, the diﬀerence between the virtual hybrid models is small: the response time history graphs virtually overlap, while the NRMSEs are all 13 DRAFT Rear and front wheel braking torqueNRMSE FOM - FFNNNRMSE FOM - PCE τrw 2.75 2.04 τfw 6.1 6.54 Table 3. NRMSE of rear and front wheel braking torque between FOM NS and RO-NSs. The values are in percent and correspondto the mean values of the 30 HS evaluations. 0 10 20 30 40051015202530 25 30 3520222426 (a)Front wheel braking torque. 0 10 20 30 4005101520 25 30 35151617 (b)Rear wheel brakingtorque. Figure 5. Front and rear wheel braking torque for the cases of FOM, FFNN-based and PCE-based MOR. Mean values of the30 HSevaluations. smallerthan2%. Therefore,boththePCEandtheFFNNMORtechniquesdescribedin Section2producevalidandeﬃcientreduced-ordermodelsforRTHS.Itisalsonotable that a combination of RO-NSsdevelopedusing diﬀerent MOR approaches also works, giving a much-neededability tocombine reduced-ordermodels inpractice. Finally, given that the performance of PCE-based and FFNN-based reduced-order models is practically the same, model designers need to make their modeling choice based on additional consideration. High-dimensional models, with long input and output vectors, are more suitable for FFNN-based MOR, since PCEs may require utilization of sparse schemes and thus incur additional errors. On the other hand, PCE-basedMORoﬀersadistinctadvantageinthatthePCEcoeﬃcientscanbeusedto compute the statistics of the relevant model response quantities, as well as assess their sensitivity to model input quantities via the Sobol’ indices approach, at no additional cost [56,57]. The latter is not a feature of neural networks. Notably, the software for developing both FFNN and PCE reduced-order models exists in the form of Matlab toolboxes [49, 33], making practical implementation of these to MOR approaches relatively easy. 5. Conclusions In this study, several model order reduction techniques are bench-marked for real- time hybrid simulation. The use of model order reduction techniques in hybrid simu- 14 DRAFT Hybrid model responsesNRMSE FOM - FFNNNRMSE FOM - PCENRMSE FOM - FFNN&PCE ωrw 1.39 1.64 1.43 ωfw 1.38 1.63 1.42 ωen 0.13 0.17 0.16 ωvd 1.39 1.64 1.43 τen 1.43 1.85 1.68 τvd 0.09 0.09 0.07 v 1.47 1.74 1.51 Table 4. NRMSE for the mean values of the selected virtual hybrid model response quantities. The values are in percent. 0 10 20 30 40-4000-3000-2000-10000 11.05 11.054-3640-3630-3620 (a)Engine crankshaft angular velocity. 0 10 20 30 40-600-400-2000200 36 38 40-100-80-60 (b)Transmission output shaft angular velocity. 0 10 20 30 40-20020406080 30.4 30.5 30.6-202 (c)Engine crankshaft torque. 0 10 20 30 40-10010203040506070 20.82 20.84 20.86-202 (d)Transmission output shaft torque. Figure 6. Selected hybrid model responses evaluated in each testing scenario. Mean values of the 30 HS evaluations. Zoomed axes are addedto highlightthe diﬀerence between the shown responsetime histories. 15 DRAFT 0 10 20 30 40050100150200250300 29 30 31120130140(a)Front wheel angular velocity. 0 10 20 30 40050100150200250300 29 30 31100120140 (b)Rear wheel angularvelocity. 0 10 20 30 4005101520253035 38 40 42234 (c)Motorcycle velocity. Figure 7. Selected hybrid model responses evaluated in each testing scenario. Mean values of the 30 HS evaluations. Zoomed axes are addedto highlightthe diﬀerence between the shown responsetime histories. 16 DRAFT lation is especially important for hybrid models encompassing high-dimensional non- linearnumericalsubstructures,riskingdistortingthetimescaleofhybridsimulationby introducingtimedelaysduetotheextendedcomputationalpowerneeded. Theselected model order reduction techniques are based on data-driven regression techniques and particularly on polynomial chaos expansion and feedforward neural networks. A para- metriccasestudyconsistingofavirtualhybridmodelisusedtovalidatetheframework. Feedforward neural networks and polynomial chaos expansion surrogates are trained to replicate the dynamic response of selected numerical substructures. Substitution of theirrespectivefull-ordercomponentsformsthereduced-orderhybridmodel. Speciﬁc parameters of the hybrid model are varied and multiple hybrid simulation evaluations are conducted to determine the robustness of the model order reduction framework. Comparisons with the full-order hybrid model responses are made to assess whether thereduced-orderhybridmodelscanaccuratelyreplicatetheinitialdynamicresponses while the corresponding errors are proven to be negligibly small. Results conﬁrm the validity, beneﬁts and performance of the proposed model order reduction framework. Inparticular,forthecasestudypresentedinthiswork,modelorderreductiontechniques allowed for a≈52%reduction of the computational cost, for a maximum of 6.54% loss on absolute accuracy versus full-order model. Model order reduction techniques are therefore capable of producing hybrid models for more reliable real-time hybrid simulations without a signiﬁcant loss of accuracy. The use of model order reduction techniques pushes the boundaries of hybrid simulation further by enabling the use of more complex and computationally involved numerical substructures. Future work will focus on quantifying how much model order reduction techniques contribute to improving real-timehybrid simulationwith physical substructures. 6. Data Availability Statement All data and code that support the ﬁndings of this study are available from the corresponding author upon reasonable request. Acknowledgments ThisprojecthasreceivedfundingfromtheEuropeanUnion’sHorizon2020research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 764547. Thesoleresponsibilityofthispublicationlieswiththeauthor(s). TheEuropean Union is not responsible for any use that may be made of the information contained herein.