Real-timehybridsimulationisanexperimentalmethodusedtoobtainthedynamic response of a system whose components consist of loading-rate-sensitive physical and numerical substructures. The coupling of these substructures is achieved by actua- tion systems, i.e., an arrangement of motors or actuators, which are responsible for continuously synchronizing the interfaces of the substructures and are commanded in closed-loop control setting. To ensure high ﬁdelity of such hybrid simulations, per- forming them in real-time is necessary. However, real-time hybrid simulation poses challengesastheinherentdynamicsoftheactuationsystemintroducetimedelays,thus modifying the dynamic response of the investigated system and hence compromising the simulation’s ﬁdelity and trust in the obtained response quantities. Therefore, a referencetrackingcontrollerisrequiredtoadequatelycompensateforsuchtimedelays. In this study, a novel tracking controller is proposed for dynamics compensation in real-time hybrid simulations. It is based on an adaptive model predictive control approach, a linear time-varying Kalman ﬁlter, and a real-time model identiﬁcation al- gorithm. Withinthelatter,auto-regressiveexogenouspolynomialmodelsareidentiﬁed in real-time to estimate the changing plant dynamics and used to update the predic- tion model of the tracking controller. A parametric virtual real-time hybrid simulation case study is used to validate the performance and robustness of the proposed control scheme. Results demonstrate the eﬀectiveness of the proposed controller for real-time DRAFT hybrid simulations. Keywords: real-time hybrid simulation, adaptive model predictive control, Kalman ﬁlter, real-time model identiﬁcation, dynamic response, actuation dynamics compensation. 1. Introduction 1 Hybrid simulation (HS), sometimes also called hardware-in-the-loop (HiL) or 2 system-in-the-loop testing, is an experimental method used to obtain the dynamic 3 response of a prototype system whose model consists of numerical (NS) and physical 4 substructures (PS). The so-called hybrid model is obtained by the coupling of NS and 5 PS.Suchcouplingisachievedusinganeﬃcientinterfacedataexchangesystemandan 6 actuationsystem,i.e.,anarrangementofmotorsoractuators,whichareresponsiblefor 7 continuously synchronizing the interfaces of the substructures and are commanded in 8 a closed-loop control setting. HS is based on a step-by-step numerical solution of the 9 equations governing the motion of the hybrid model [1]. In more detail, in each time 10 stepofHS,theNSgeneratesacommandthatneedstobetrackedbythePStopreserve 11 synchronization of the substructures’ boundary conditions. The dynamic response of 12 the PS is then measured using appropriate sensors and fed back to the NS to establish 13 the state of the hybrid modeland advancethe HSto thenexttime step. 14 HSisadesirabledynamicresponsetestingtechniquesinceitcombinestheversatility 15 and risk-free testing of numerical simulations along with the realism of experimental 16 campaigns. However, ensuring high ﬁdelity of hybrid dynamic response simulations 17 of a loading-rate-sensitive prototype often necessitates performing HS in real-time. 18 Real-time hybrid simulation (RTHS) poses several challenges. One challenge arises 19 from complex NS because of the computational power required to compute the NS 20 response on time, i.e., within the HS time step. For example, to adequately capture 21 the dynamic response of high dimensional and/or nonlinear numerical models, it is 22 common practice to decrease the time step of the NS simulation which increases the 23 demandedcomputationalpower. Notcompletingthecomputationstoestablishthestate 24 of the NS on time introduces delays and risks distorting the time scale of HS. In such 25 2 DRAFT cases, model order reduction of the NS is a very eﬀective countermeasure [2, 3]. 26 Another challenge in RTHS originates from the PS, speciﬁcally from the actuation 27 system,asitsinherentdynamicsintroducetimedelaysthatmodifythedynamicresponse 28 ofthePSandhencecompromisetheﬁdelityoftheentireHS.Addressingthischallenge 29 is the focus of this paper. 30 VariouscontrolschemesfordynamicscompensationinRTHSexistintheliterature. 31 Aselectionofsuchschemesispresentedbelow. Aninitialapproachhasbeendeveloped 32 by Horiuchi using a polynomial extrapolation algorithm [4], which was afterwards 33 extended into an adaptive scheme [5, 6]. Compensation for the phase shift of the 34 actuation system, via phase-lead compensators, was also proposed by various authors 35 [7,8,9]. InversecompensationisanothertechniqueusedinRTHS,inwhichaninverse 36 model of the actuation system is employed as a feedforward compensator [10, 11]. In 37 ordertoimprovetherobustnessofRTHSinthepresenceofexperimentaluncertainties, 38 adaptive compensation strategies were developed aiming at estimating in real-time the 39 controller parameters [12, 13,14,15]. 1loop shaping control schemes were also 40 proposed as an alternative robust approach [16, 17]. Additionally, a feedforward- 41 feedback control scheme based on linear-quadratic regulator (LQR) and Kalman ﬁlter 42 wasdevelopedforcontrollingbothsingleandmulti-actuationsetups[18, 19]. Amodel- 43 based sliding mode control approach has also been developed for RTHS making use 44 of a reduced plant [20]. Recently, a control scheme based on model predictive control 45 (MPC) [21]was alsopresented. 46 Asshownin[21],MPCisaneﬀectivecontrolstrategyforRTHS.Akeyadvantageof 47 MPCisitsabilitytoperformoptimizationinreal-timeaswellassatisfyingconstraintsof 48 theexaminedsystem,e.g.,actuationsystemcapacitylimitations. However,inMPCthe 49 real-timeoptimizationisperformedbasedonpredictionsofthefuturedynamicresponse 50 of the system under control. Therefore, the performance of the optimization depends 51 on the accuracy of the predictions. Classical MPC utilizes a prediction model, e.g., a 52 lineartime-invariant(LTI)dynamicmodel,forthesepredictions. Incasethedynamics 53 of the examined system are highly nonlinear or linear with time-varying parameters, 54 the predictions from a LTI model may not be accurate and thus the performance of 55 MPC will be degraded. In this regard, an adaptive MPC is a suitable solution, since 56 3 DRAFT theLTImodelusedtopredictthedynamicresponseofthecontrolledsystemisadapted 57 during the simulation in order to capture the changing dynamics of the system under 58 control (the co-called plant). Some of the common practices for model updating in 59 MPC include real-time systemidentiﬁcation andsuccessive linearization [22, 23]. 60 ThisstudypresentsanoveltrackingcontrollerfordynamicscompensationinRTHS 61 extending the classical MPC of [21] to an adaptive scheme. It is based on the adaptive 62 MPCapproach,alineartime-varyingKalmanﬁlter,andareal-timemodelidentiﬁcation 63 algorithm. Within the latter, auto-regressive exogenous (ARX) polynomial models are 64 identiﬁed in real-time to estimate the changing plant dynamics and used to update the 65 prediction model of the adaptive MPC. A parametric virtual RTHS case study is used 66 todemonstrateandvalidatetheperformanceandrobustnessoftheproposedcontroller. 67 The paper is organized as follows. Section 2describes the tracking controller. 68 Section3introduces the virtual RTHS case study where the proposed controller is 69 developed, validated and demonstrated. Section 4discusses the obtained results and 70 Section5presentsthe overall conclusionsof thisstudy. 71 2. The tracking controller 72 In this section, the individual components of the proposed tracking controller are 73 described. Asmentionedabove,thesecorrespondtotheadaptiveMPC,thelineartime- 74 varyingKalmanﬁlterandareal-timemodelidentiﬁcationalgorithmusedtoupdatethe 75 prediction model of the MPC. Figure 1illustrates the block diagram of the controller. 76 The plant corresponds to the system under control. The latter is addressed in more 77 detail inSection 3. 78 2.1. Adaptivemodel predictive control 79 In the MPC approach, the current control action is calculated by solving, in each 80 control interval, a ﬁnite horizon optimization problem utilizing the current state of the 81 plant as its initial conditions. The outcome of this optimization is an optimal control 82 sequence. Theﬁrstcontrolactionofthissequenceisusedinthecurrentcontrolinterval 83 [24,25]. Controlintervalsaresetsofconsecutivetimesteps,obtainedfromasampling 84 4 DRAFT Tracking ControllerAdaptive MPCPlant Kalman Filtery Measured System Outputr y^zn Real-TIme Model Identification zyReferenceDisturbance dNoiseFigure 1. Architecture of the proposedtracking controller. frequency that the MPC uses internally. Essentially, a MPC solves a ﬁnite horizon 85 optimization problem similar to 1and LQR control techniques. A fundamental 86 diﬀerence nevertheless lies in the fact that MPC solves it in real-time; namely, at each 87 control interval a new optimization problem in a receding horizon approach is solved. 88 Therefore,thecontrolactionsareupdatedinreal-time,whileinclassicaloptimalcontrol, 89 e.g., LQR, the optimization is solved oﬄine and thus the control law is deﬁned prior 90 to the simulation and remains constant. Additionally, MPC oﬀers the capability of 91 satisfyinghardconstraintsregardingtheinputs,outputsand/orstatesoftheplantwhile 92 solving the optimization problem in real-time [26]. 93 AMPCencompassesfourelements,namelythepredictionmodel,theperformance 94 index (or the cost function), the constraints, and a solver used to compute the control 95 actions. The prediction model constitutes the core of the MPC as it is used to predict 96 the future plant outputs based on information available up to the current time instant, 97 e.g., previous and current inputs/outputs, and on the future control actions [27]. As 98 mentionedabove,theoverallperformanceofMPCisdirectlyinﬂuencedbytheresponses 99 ofthepredictionmodel. Therefore,thismodelshouldcapturethedynamicsoftheplant 100 as accurately as possible. Although the MPC updates the control actions in real-time, 101 the controller is still not adaptive in its structure, e.g., the prediction model remains 102 the same. Updating the prediction model in real-time results in the so-called adaptive 103 MPC. 104 5 DRAFT Plant Prediction Model Optimizer Cost functionConstraintsPredictionsControl sequence ReferenceMPC ControllerOutputs Inputs y^ rzk J ∗(zk) Ax ≤ bFigure 2. MPC structure ofthe trackingcontroller. Figure 3. MPC methodology. Figure2displaysthestructureoftheMPCandFigure3 illustratesitsmethodology. 105 In more details, the prediction model predicts at each control interval :the next plant 106 outputs ˆH¹:¸8j:º 8=0%for every time instant 8within a prediction horizon %. 107 The notation ˆH¹:¸8j:ºcorrespond to the output at the time instant :¸8evaluated at 108 :. The control sequence I) :=»D¹:j:º)D¹:¸8j:º)D¹:¸%1j:º)¼includes 109 thesequenceoffuturecontrolactions D¹:¸8j:º. Ateach:,anewcontrolsequenceis 110 computed optimizinga newcostfunction expressedby: 111 J¹A:ˆH:I:º==HÕ 9=1%Õ 8=1 FH9»A9¹:¸8j:ºˆH9¹:¸8j:º¼ 2¸ (1) =DÕ 9=1%1Õ 8=0 FD9»D9¹:¸8j:ºD9¹:¸81j:º¼ 2 where=Histhenumberofplantoutputs, =Dthenumberofplantinputs, A9¹:¸8j:ºthe 112 reference value to be tracked at the 8-th time instant step from the 9-th plant output, 113 ˆH9¹:¸8j:ºthepredictedvalueofthe 9-thplantoutputatthe 8-thtimeinstant, D9¹:¸8j:º 114 the9-thplantinputatthe 8-thtimeinstant, FH9thetuningweightofthe 9-thplantoutput 115 andFD9the tuning weight of the 9-th plant input. The %,FH9andFD9parameters of 116 Eq. (1)are computedoﬄine andremain constant inthe workpresented herein. 117 In addition, an output disturbance 3and a measurement noise model =are used in 118 theproposedtrackingcontroller. Theyareimplementedasadditivetotheplantoutputs, 119 representingpotentialdisturbancesandsensornoisethatcouldbepresentinRTHS.The 120 6 DRAFT models of the disturbance 3and the noise =are expressed in state-space, as described 121 in Eqs. (2) and (3), respectively: 122 G3¹:¸1º=3G3¹:º¸3D3¹:º (2) 3¹:º=3G3¹:º¸3D3¹:º G=¹:¸1º==G=¹:º¸=D=¹:º (3) =¹:º==G=¹:º¸=D=¹:º where3,3,3,3and=,=,=,=arematricesassociatedwiththedisturbance 123 3and noise=respectively. The inputs to these state-space models, namely D3for the 124 disturbance and D=for the noise, are Gaussian variables with zero mean and unit 125 variance. 126 The control sequence I:is computed at the beginning of each control interval in 127 theoptimizer(seeFigure 2). ThecostfunctionofEq.(1)isformulatedintoaQuadratic 128 Programming (QP) problem [28, 29] and this is what the optimizer eventually solves. 129 The QPformulationadmits: 130 min G¹1 2G) G¸5)Gº (4) subject to G1 (5) where theG1inequality represents the applied constraints, Gthe solution vector 131 and itsHessianmatrix. Thevectors and1deﬁnetheconstraints,and 5isavector 132 computed using: 133 5= G¹:j:º)G¹:j:º¸ A¹:j:º)A¹:j:º¸ D¹:j:1º)D¹:j:1º(6) whereA¹:j:ºis the reference signal at the current control interval, D¹:j:1ºis the 134 applied control action in the previous control interval and a weighting factor. For a 135 more comprehensive review ofQP,the readeris encouraged to consult [30, 31]. 136 7 DRAFT In the proposed tracking controller, the QP problem of Eqs. (4) and (5) is solved 137 employing an active-set solver which applies the KWIK algorithm [32]. The latter is a 138 built-in QP solver from the Model Predictive Toolbox of MATLAB, used in this study 139 to compute the control law sequence [33]. 140 The utilized adaptive MPC algorithm is summarized as follows: 141 1.Prediction step: Considering a discrete time multiple-input-multiple-output 142 (MIMO) system that represents a linearized model of the plant, the prediction 143 model in state-space formulation at the control interval :follows: 144 G?¹:¸1º=: ?G?¹:º¸: ?D¹:º (7) H?¹:º=: ?G?¹:º¸: ?D¹:º where: ?,: ?,: ?and: ?are matrices associated with the prediction model at 145 the control interval :, representing the plant dynamics. The disturbance and 146 noise models of Eqs. (2) and (3) respectively, are assumed to be additive to 147 the measured plant output Hsignal in Figure 1, and hence the latter admits: 148 H¹:º=H?¹:º¸3¹:º¸=¹:º. 149 2.Identiﬁcation step: At each control interval :, an ARX polynomial model of 150 the plant is identiﬁed and used to adapt the prediction model of MPC, i.e., the 151 matrices: ?,: ?,: ?and: ?, (see Section 2.3for a more detail description of 152 the algorithm). 153 3.Updatestep: Ateachcontrolinterval :,the!and"vectorsoftheKalmanﬁlter 154 areupdatedtobeconsistentwiththeadaptedpredictionmodel(seeSection 2.2). 155 4.Optimizationstep: Usingthecurrentstateestimatesandpredictedplantoutputs 156 from the Kalman ﬁlter (Eqs. (16) and (19)), the MPC solves the optimization 157 problem described byEq. (1)for thecurrentcontrol interval :: 158 min I:J¹A:ˆH:I:º (8) subject to D9<8=D9D9<0G (9) 8 DRAFT where the above constraints represent safety limits regarding input signals to the 159 actuation system. 160 5.Controlstep: Thecontrolsequence I) :=»D¹:j:º)D¹:¸8j:º)D¹:¸% 161 1j:º)¼for the current control interval is computed and its ﬁrst control action is 162 applied to the plant. 163 6. Steps 1-5 are repeated till completion of the RTHS. 164 If the control interval coincides with the sampling frequency of RTHS, then at 165 the next time step of RTHS a new optimization will be performed and a new control 166 sequence will be computed taking into account the new information available, e.g., 167 disturbance,noise,introducedinthattimestep. Whenthecontrolintervalislargerthan 168 the RTHS time step, i.e., includes several time steps, then the ﬁrst control action of 169 step5willbeappliedaslongastheongoingcontrolintervallasts. Therefore,reducing 170 the time steps in the control interval would result in a more robust MPC. However, the 171 smallerthecontrolintervalthemoretheoptimizationsthatmustbeperformed,thusin- 172 creasingthecomputationalload. Inthisregard,ahigherdimensionalpredictionmodel, 173 that could potentially capture more accurately the plant dynamics, would increase the 174 performanceofMPCasthepredictionswouldbemoreaccurate. Nevertheless,ahigher 175 dimensionalmodelwouldrequireadditionalcomputationalpower. Consequently,taken 176 the aforementioned points into consideration, there existatrade-oﬀ between controller 177 accuracy and its computational performance and simulation model accuracy and its 178 computational performance. 179 2.2. Linear time-varyingKalman ﬁlter 180 As mentioned above, the performance of the MPC depends on the accuracy of the 181 predicted plant outputs. In order to provide the MPC with low-noise plant outputs, a 182 Kalman ﬁlter is used. Furthermore, the Kalman ﬁlter provides estimates of the plant 183 outputswhentherearenoneavailablefromtheplantsensors. Sincethepredictionmodel 184 of the MPC is adapted in real-time, the gain vectors of the Kalman ﬁlter also need to 185 be adapted in real-time, resulting in a linear time-varying Kalman ﬁlter. Considering 186 Eqs. (7), (2) and (3), the state-space formulation including the plant dynamics, the 187 disturbance and noise model follows: 188 9 DRAFT G2¹:¸1º=: 2G2¹:º¸: 2D2¹:º (10) ˆH¹:º=: 2G2¹:º¸: 2D2¹:º where 189 : 2=266666664: ?0 0 030 0 0=377777775 : 2=266666664: ?0 0 030 0 0=377777775 (11) : 2=h : ?3=i  : 2=h : ?3=i  G) 2=h G) ?G) 3G) =i andD) 2=h D)D) 3D) =i  : ?,: ?,: ?and: ?matrices are updated in real-time while 3,3,3and3and 190 =,=,=and=of Eqs. (2) and (3) are considered constant during RTHS. The G2, 191 D2vectors contain the states and inputs of the plant ?, the disturbance 3and noise 192 =, respectively. The weighting coeﬃcients &'#of the Kalman ﬁlter are constant 193 during RTHS and are computed from the following expectations: 194 &=E 33)  ' =E ==)  # =E 3=)  (12) Figure1displaystheinterconnectionoftheKalmanﬁlterwithintheproposedtracking 195 controller. In more detail, at each control interval :, theG2states are estimated as 196 follows: 197 1.Gain computations: The gain vectors !"of the Kalman ﬁlter are updated to 198 be consistent withthe adaptedprediction model: 199 !:= : 2%:j:1:) 2¸# : 2%:j:1:) 2¸'1  (13) ":=%:j:1:) 2 : 2%:j:1:) 2¸'1  (14) %:¸1j:=: 2%:j:1:) 2 : 2%:j:1:) 2¸# !:)¸& (15) 10 DRAFT where!:":denote the!"gain vectors at the control interval :and%:¸1j: 200 the state estimate error covariance matrix at :¸1, calculated with information 201 obtained fromthe :control interval. 202 2.Measurement update step: The current state estimate G2¹:j:ºis adjusted with 203 the latest measurements: 204 G2¹:j:º=G2¹:j:1º¸": H¹:º: 2G2¹:j:1º  (16) whereG2¹:j:1ºdenotes the state estimate from the :1control interval. 205 3.Prediction step: The state for the next, :¸1, controlintervalis estimatedas: 206 G2¹:¸1j:º=: 2G2¹:j:1º¸: 21D¹:º¸!: H¹:º: 2G2¹:j:1º (17) whereD¹:ºthe control action used from ¹:1ºtill:,H¹:ºthe plant output 207 measured at :and: 21the ﬁrst columnof : 2. 208 Oncethecurrentstateoftheplantisestimated,thecontrollerpredictsplantoutputsfor 209 the entireprediction horizon ofthe current controlinterval. 210 1. For any successive time instant of the prediction horizon, 8=1%, the 211 predicted plantstatesare obtainedby: 212 G2¹:¸8j:º=: 2G2¹:¸81j:º¸: 21D¹:¸81j:º (18) 2. The predictedplant outputfor 8=1%is obtainedby: 213 ˆH¹:¸8j:º=: 2G2¹:¸8j:º (19) 2.3. Real-timemodel identiﬁcationalgorithm 214 At each control interval :, a real-time model identiﬁcation algorithm is employed 215 and ARX polynomial models [34] of the plant are identiﬁed and then used to update 216 the prediction model of the MPC. The formulation of ARX polynomial models is the 217 following: 218 11 DRAFT ¹@ºH¹Cº=¹@ºD¹Cº¸4¹Cº (20) where@is a time shift operator, ¹@º=1¸U1@1¸U2@2¸¸U=0@=0,¹@º= 219 11¸12@1¸13@2¸¸1=1@¹=11ºandD¹Cº,H¹Cºand4¹Cºtheinputs,outputsand 220 error respectively. 221 The model identiﬁcation is based on the recursive inﬁnite-history algorithm [35] 222 following: 223 \¹Cº=\¹C1º¸ ¹CºH¹CºHpred¹Cj\º (21) where 224 \¹Cº=h U1¹CºU2¹Cº U=0¹Cº11¹Cº 11=¹Cºi (22) corresponds to the identiﬁed parameters at time C,H¹Cºthe measured plant output at C, 225 Hpred¹Cj\º=k)¹Cº\¹C1ºthe prediction of H¹Cºaccounting for measurements up to 226 C1and ¹Cº=&¹Cºk¹Cº, wherek¹Cºis thegradient of Hpred¹Cj\ºand&¹Cºadmits: 227 &¹Cº=%¹C1º 1¸k)¹Cº%¹C1ºk¹Cº (23) where 228 %¹Cº=%¹C1º¸1%¹C1ºk¹Cºk)¹Cº%¹C1º 1¸k)¹Cº%¹C1ºk¹Cº (24) %¹Cºrepresents the covariance ofparameter identiﬁcationerrorwith %¹0º=1. 229 At each control interval, the parameters \¹Cºare identiﬁed and the ¹@º,¹@º 230 polynomialareconvertedtostate-spaceformulationresultinginthe : ?,: ?,: ?and: ?231 matrices used in Eq. (7). 232 Asmentionedearlier,aKalmanﬁlterisusedinthecontrolschemeforstateestima- 233 tionandtherecursiveinﬁnite-historyalgorithmforreal-timeparameteridentiﬁcationof 234 the prediction model. As displayed in Figure 1, these are implemented as two separate 235 blocks, in a decentralized approach. This approach oﬀers the ﬂexibility of the inde- 236 pendent design of the state and the parameter estimation algorithms. Consequently, a 237 12 DRAFT varietyofcombinationsofdiﬀerentreal-timemodelidentiﬁcationalgorithmscaneasily 238 be tested. A single joint algorithm, a centralized approach, featuring both state and 239 parameterestimationcouldalsobeimplemented,butthisisnotexaminedinthisstudy. 240 3. Case study 241 To validate the performance and robustness of the proposed AMPC tracking con- 242 troller, a parametric case study is utilized. To facilitate the design and testing of the 243 controller,thehybridmodelissimulatedvirtuallyinthesensethatallofitssubstructures 244 are simulated numerically and thus virtual PS (vPS) are employed instead of physical 245 specimens. In this section, the utilized performance metrics are deﬁned ﬁrst, then the 246 problem formulation of the case study is introduced and, ﬁnally, the design properties 247 of the tracking controller are addressed. 248 3.1. Performance metrics 249 To assess the performanceof the controller,three metricsare deﬁned. Theseare: 250 1. Tracking time-delay, deﬁned by the correlation between the reference signal A 251 and plant output H: 252 1= arg max :Corr¹A¹8ºH¹8:ºº 5RTHS[msec] (25) where5RTHS=10kHz is the sampling frequency of RTHS. 253 2. NormalizedRootMeanSquareError(NRMSE)ofthetrackingerror,denotedas: 254 2=vuuuuutÍ# 8=1h H¹8ºA¹8ºi2 Í# 8=1h A¹8ºi2100[%]. (26) 3. Peak TrackingError(PTE), denotedas: 255 3=maxjH¹8ºA¹8ºj maxjA¹8ºj100[%]. (27) 13 DRAFT Above, 1represents the maximum cross-correlation between the reference and the 256 measured plant output. It is used to quantify how these two signals diﬀer in time. 257 Speciﬁcally, cross-correlation represents the number of time steps that the measured 258 plant output should be shifted to match the reference. For 1¡0the measured plant 259 output has a delay with respect to the reference, i.e., there is a tracking delay, while 260 for 10themeasuredplantoutputisleadingthereference,i.e.,overcompensationis 261 occurring. Performancemeasure 2establishesthelevelofhowquantitativelydiﬀerent 262 the reference and measured plant output are, considering the entire simulation period. 263 Performance measure 3represents the maximum tracking error. The objective of the 264 tracking controller is to maintain 1, 2and 3as close to zero as possible [36]. 265 3.2. Problem formulation 266 The case study prototype is a motorcycle. The hybrid model of the motorcycle 267 is made of four NS and one vPS. The NS are: i) the engine, ii) the motorcycle body 268 dynamics,iii)therearwheelbrakingsystemandiv)thefrontwheelbrakingsystem. The 269 vPS is the electrically continuouslyvariable transmission(eCVT) ofthe motorcycle. 270 The eCVT vPS corresponds to a MIMO model with two sets of one input / one 271 output. TheﬁrstsetisconnectedtothemotorcycleengineNSandthesecondonetothe 272 motorcyclebodyNS.Thelatterconnectioncorrespondstothetransmissionoutputshaft 273 ofthemotorcycle. TheengineNSsimulatesthedynamicsofthecombustionengine. It 274 is represented by a multi-input-single-output (MISO) model, with its inputs being the 275 throttlepercentage CAandtheangularvelocityoftheengine l4=anditssingleoutput 276 being the torque of the engine g4=. The motorcycle body NS represents the body or 277 chassisdynamicsofthemotorcycletogetherwiththedynamicsofthewheels,tiresand 278 suspensionsalongwiththeroadproﬁleandtheenvironmentaldrivingconditions. Itis 279 represented by a MIMO model with 3 sets of one input / one output. The ﬁrst set is 280 connected to the eCVT vPS with the torque gE3as input and the angular velocity lE3281 of the transmission output shaft as output. The second and third sets are linked to the 282 rear and front wheel braking system NS respectively. Both braking systems are MISO 283 models. The rear wheel braking system NS inputs are the angular velocity of the rear 284 wheellAFandtheappliedforceonthebrakepedal 1AAFwhileitssingleoutputisthe 285 14 DRAFT τfw ωfwτrwτenref ωrwωvdωenNSEngineeCVT vPS Motorcycle body dynamics NS Front wheel braking system NS Rear wheel braking system NSτvdTracking Controllerτen Fbrfwthr FbrrwGmInverter + MotorFigure 4. Hybrid modelblock diagram. brakingtorqueoftherearwheel gAF. Respectively,thefrontwheelbrakingsystemNS 286 inputsaretheangularvelocityofthefrontwheel l5 Fandtheappliedforceonthebrake 287 lever1A5 Fanditssingleoutputisthebrakingtorqueinthefrontwheel g5 F. Figure4 288 illustratestheinterconnectionsbetweenNSandthevPS,whileFigure 5depictsthereal 289 eCVTPSthatwouldbeutilizedinanon-virtualRTHS.Todevelopthesubstructuresof 290 the hybrid model, the multi-physics simulation software Simcenter Amesim was used. 291 Thereportof[37]oﬀersathoroughdescriptionoftheutilizedsubstructuresalongwith 292 theequationsgoverningtheirmotion. Tointerconnectandco-simulateallsubstructures 293 aswellastheRTHSalgorithminreal-time,the Simcenterreal-timeplatform wasused. 294 The forces applied on the brake pedal (to activate the rear wheel brakes) and the 295 brake lever (to activate the front wheel brakes) are considered equal and expressed in 296 Newton,1AAF=1A5 F, andCA,1AAFas deﬁned in Eqs. (28) and (29) respectively. 297 Themaximumappliedthrottleis0.5(50%),whichcorrespondstohalf-openthrottle. In 298 Figure6the applied driving scenario of the case study is shown. The road proﬁle, i.e., 299 height of the ground, is expressed by the ¹Gºsinusoidal signal and follows Eq. (30), 300 whereGis the current motorcycle position in meters. The ambient wind velocity is 301 considered to be zero. The simulation duration of the case study is 45 sec. A fourth- 302 order Runge–Kutta (RK4) method with a ﬁxed time step of 0.1 msec is used as the 303 numerical integration scheme in the conducted RTHS of the dynamic response of the 304 motorcycle virtual hybrid model during the driving scenario. Figure 7displays some 305 indicative dynamic responses of the virtual hybrid model of the motorcycle, excited 306 15 DRAFT Figure 5. The eCVT testbench. 0 10 20 30 4000.10.20.30.40.5 01020304050 Figure 6. Case study drivingscenario. under the driving scenario on the given roadproﬁle andin thegiven windconditions. 307 CA¹Cº=8>>>>>>>>>> < >>>>>>>>>>:0125C,5C 9 05,9C13 0125C,13C17 0 , elsewhere(28) 1AAF¹Cº=8>>>>>>>>>> < >>>>>>>>>>:10C,20C 25 50,25C32 10C,32C37 0, elsewhere(29) ¹Gº=8>>> < >>>:0 ,0G2 002 cos2cG 31 , elsewhere(30) 16 DRAFT 0 10 20 30 40050100150200250300(a) 0 10 20 30 4005101520 (b) 0 10 20 30 40-600-400-2000200 (c) 0 10 20 30 40-100102030405060 (d) 0 10 20 30 40-505101520 (e) 0 10 20 30 4005101520253035 (f) Figure 7. Indicative hybrid model responses: (a)angular velocity of the rear wheel, (b)braking torque of the rear wheel, (c)angular velocity of the transmission output shaft, (d)torque of the transmission output shaft,(e)torque of theengine and (f)motorcycle velocity. BecausetheeCVTvPS(themotorcycletransmission)hastwoinputs, g4=andlE3, 308 two actuation systems would be required if the RTHS would be conducted physically; 309 17 DRAFT hence two individual tracking controllers would be needed. However, without loss 310 of generality and for the sake of simplicity, one actuation system/tracking controller is 311 utilizedasillustratedinFigure 4. ItisappliedbetweentheeCVTvPSandtheEngineNS 312 astheinterconnectionbetweenthetwohasthehighestfrequencycontentthusrequiring 313 the compensationfor theRTHS tobe asgood aspossible. 314 TheplantofthecasestudycorrespondstotheeCVTvPSactuationsystem,i.e.,the 315 inverter and the electric motor of Figure 4. Therefore, the plant output Hof Figure 1 316 corresponds to the torque of the engine g4=while the reference Asignal tog4=ref. 317 Measurements were conducted on the corresponding test bench at the testing facilities 318 ofSiemensIndustrySoftwareforthesystemidentiﬁcationoftheinverterandtheelectric 319 motor. This actuation system was identiﬁed as a second-order transfer function. Such 320 system isrepresented incontinuous timeas follows: 321 <=11B¸10 02B2¸01B¸00 (31) where 322 11=0 1 0=00=5657 0 2=1and01=7862 (32) It can also be expressed as <=l2 =¹B2¸2Zl=B¸l2 =º, withl=andZits natural 323 frequency and damping, respectively. 324 In a study by Silva et al. [36], coeﬃcients of variation of up to 5% and 7% of the 325 plant’s natural frequency and damping, respectively, were identiﬁed after conducting 326 severalHSruns. Thesewere mainlyduetophysicalchanges oftheexperimentalsetup 327 after repeated testing. Therefore, in order to simulate potential variations of the plant 328 <usedinthisstudy,itsnaturalfrequency l=anddamping Zwereselectedtovary. To 329 beconservativeandexposethetrackingcontrollertomoresevereconditions, l=andZ 330 varywithacoeﬃcientofvariationof15%and17%respectively,duringthesimulation, 331 every 0.25 seconds of the RTHS. The variations are implemented by changing the 332 valuesof10,00. Particularly, 10,00N¹ 5657 164053ºand10equals00. Therefore, 333 the goal of the real-time model identiﬁcation algorithm, presented in Section 2.3is to 334 capture the applied variations of <and to correctly identify its changing dynamics 335 18 DRAFT andprovidetheMPCwiththeupdatedpredictionmodel. Figures 9a-dillustrateinblue 336 color the applied variations of <over the simulation time. The initial values of <337 are thoseof Eq.(32). 338 Furthermore, to assess the robustness of the proposed tracking controller under 339 varioushybridmodelconﬁgurations,elevenparametersofthemotorcyclevirtualhybrid 340 modelareselectedtovary. Table 1providesanoverviewoftheseparameters. Auniform 341 distribution is assigned to each chosen parameter. Their mean values and standard 342 deviations are identiﬁed from [38, 39,40] to reﬂect a range of possible parameter 343 variationsofthecorrespondingmotorcyclecomponents. Thenominalparametervalues 344 for this case study are their mean values. Using the Latin hypercube sampling (LHS) 345 method [41], 200 samples are obtained from each parameter and hence 200 hybrid 346 model evaluations are performed. In contrast with 0010parameters of <that are 347 varyinginreal-time,theparametersofTable 1remainconstantduringeachRTHSrun, 348 but their values are changed in every successive RTHS. For the scope of the presented 349 study,200samplepointsareconsideredsuﬃcienttoadequatelycapturetheunderlying 350 probabilityspaceandthustoexposethetrackingcontrollertowideparametervariations. 351 Thestabilityofthehybridmodelwasnotaﬀectedbythevariationoftheplantdynamics 352 andTable’s 1parameters. ThiscanbeappreciatedfromtheresultspresentedinTable 2 353 and Figure 8of Section 4, asin each casethe hybrid modelresponse was stable. 354 3.3. Tracking controller designproperties 355 Theinitialmatricesofthepredictionmodel( 0 ?0 ?0 ?and0 ?ofEq.(7))coincide 356 with thedynamics ofthe plantmodel ofEq. (31)and follow: 357 0 ?=26666409685900348 002519 099963777750 ?=26666400031 0377775 (33) 0 ?=h 0 110488i 0 ?=h 0i  The disturbance 3and noise=models, described in Eqs. (2) and (3), respectively, are 358 additive tothe plantoutput andadmit: 359 19 DRAFT Parameter Description DistributionMean valueStandard deviationUnits ACVertical stiﬀness rear tire Uniform 58570 11714 N/m /ACVertical damping rear tire Uniform 11650 3495 Ns/m 5CVertical stiﬀness front tire Uniform 25000 5000 N/m /5CVertical damping front tire Uniform 2134 6402 Ns/m ABStiﬀness rear suspension Uniform 125000 25000 N/m /ABDamping rear suspension Uniform 10000 3000 Ns/m 5BStiﬀness front suspension Uniform 19000 3800 N/m /5BDamping front suspension Uniform 1250 375 Ns/m "Motorcycle massUniform 300 6 Kg Engine moment of inertia Uniform 00115 00023 Kgm2 `Engine coeﬃcient of viscous friction Uniform 0001 000005 Nm/(rev/min) Table 1.Motorcycle virtual hybrid modelparameters withtheirassigned statistical distributions. 20 DRAFT G3¹:¸1º=G3¹:º¸00004D3¹:º (34) 3¹:º=G3¹:º =¹:º=D=¹:º (35) with3=13=000043=1and3=0, while==0==0== 360 0and==1. Recall that the inputs to the disturbance and noise models, i.e., D3and 361 D=, respectively, follow the standard Gaussian distribution, namely D3D=N¹ 01º. 362 According to Eqs. (10) and (11), the state-space formulation of the plant dynamics for 363 :=0, includingthe disturbance andnoise modelfollows: 364 2=26666666666409685900348 0 0 002519 09996 0 0 0 0 1 0 0 0 0 03777777777752=26666666666400031 0 0 0 0 0 0 00004 0 0 0 0377777777775(36) 2=h 0 110488 1 0i 2=h 0 0 1i  The initialKalman ﬁltergainvectors admit: 365 !0=h 0157 0917 0398i) 1e3 (37) "0=h 0195 0912 0398i) 1e3 The tuning weights of Eq. (1) are selected to be FH=5204andFD=0096, while 366 =D=Hequalto1astheplantisasingle-input-single-output(SISO)model. Eachcontrol 367 intervalisobtainedwithasamplingfrequencyof2.5kHzandthepredictionhorizonis 368 %=10. The constraints of Eq. (9) represent safety limits regarding the plant input and 369 admit: 370 200D200[Nm] (38) 21 DRAFT 4. Results and discussions 371 Table 2presents the performance metrics, 1, 2and 3, of the proposed tracking 372 controller for the aforementioned simulation and Figure 8displays their normalized 373 histograms. 374 NominalStochastic Mean Values Standard Deviation 1[msec] 2.5 2.5 0 2[%] 0.32 0.34 0.01 3[%] 2.79 2.82 0.21 Table 2.Performance metrics results. Figure 8. Normalized histograms of 1, 2and 3, obtained from200 virtualRTHSevaluations. ThenominalresultsofTable 2refertothecasewhenthemeanvaluesoftheparameters 375 in Table 1were used: hence they correspond to a single deterministic RTHS. The 376 stochastic results refer to the outcomes from the 200 RTHS evaluations. In both 377 deterministic (nominal) and stochastic cases, the parameters 00,10of <vary as 378 displayed inFigure 9. 379 As mentioned above, the performance of the tracking controller is assessed by 380 how close to zero the performance metrics are. From Table 2and Figure 8, it can be 381 22 DRAFT 0 10 20 30 4002000400060008000100001200014000(a) 0 10 20 30 40747678808284 (b) 0 10 20 30 40020004000600080001000012000 (c) 0 10 20 30 40-0.500.511.522.5 (d) Figure9. Appliedvariationsintheparametersoftheplantandtheidentiﬁedvaluesfromthereal-timemodel identiﬁcation algorithm: (a)00,(b)01,(c)10and(d)11. appreciatedthat 1, 2and 3arequitesmallandhencetheproposedcontrolschemecan 382 adequately regulate the desired plant output even under the presence of hybrid model 383 parameter and plant dynamics variations,proving aswell itsrobustness. 384 Figure 9 also presents in red color the parameters of the plant dynamics that the 385 real-timemodelidentiﬁcationalgorithmestimatesateachcontrolinterval. Asstatedin 386 Eq. (28), for the ﬁrst 5 seconds of the simulation there is no input (zero values) to the 387 hybridmodelandhencetheoutputisalsozero. Asaresult,theidentiﬁcationalgorithm 388 returnszerovalues. Itcanbeobservedthataftertheﬁrst5seconds,thealgorithmstarts 389 to respond. Recall that there were no applied variations for 01and11. However, the 390 algorithmfailstoaccuratelyidentifythesetwoparameters(Figures 9b,d). Nevertheless, 391 theirdeviationisrelativelysmall. Regarding 00and10,thepredictionsaresatisfactory 392 23 DRAFT 0 10 20 30 40-0.02-0.0100.010.020.03Figure 10. Plant system identiﬁcation error. 0 10 20 30 40-505101520 11.88 11.9 11.9218.218.418.6 20.1 20.2 20.3-0.6-0.4-0.2 Figure 11. Torque responses ofthe engineNS: referenceg4=A4 5and measured g4=. as shown from Figures 9a,c. Figure 10 displays the overall identiﬁcation error of the 393 real-time model identiﬁcation algorithm, deﬁned as: 394 Y¹Cº=H¹CºHpred¹Cº maxjH¹Cºj100»%¼ (39) FromFigure10,itcanbestatedthattheoverallidentiﬁcationerrorisquitesmallandthus 395 the parameter identiﬁcations are acceptable. It should be noted that the performance 396 ofthereal-timemodelidentiﬁcationalgorithmiscrucialfortheoverallperformanceof 397 theadaptiveMPC.ThemodelthatthealgorithmidentiﬁesisusedinternallyinMPCas 398 its prediction model. Hence, faulty system identiﬁcation of the plant could yield large 399 tracking errors. 400 Figure 11 displays, for the nominal RTHS case, the reference values of the plant, 401 namely the output of the engine NS g4=A4 5, and the plant output, namely the input 402 to the vPS, g4=. An ideal tracking controller should be able to compensate for the 403 actuation system dynamics so that g4=A4 5andg4=are as close as possible. As shown 404 from Figure 11, the tracking error betweenthe twosignals isverysmall indeed. 405 Basedontheaboveresults,theproposedtrackingcontrollercanprovidethedesired 406 performance and at the same time be robust to hybrid model and actuation system 407 variations,aswellastotheintroduceddisturbancesandmeasurementnoise. Addition- 408 ally,usingtheproposedtrackingcontroller,thetrackingerrorsandtimedelaysthatare 409 introducedduetotheinherentdynamicsoftheactuationsystemused,aresatisfactorily 410 24 DRAFT compensated,enablingthusRTHSoutcomesofhighﬁdelity. Asaresult,theproposed 411 control scheme is demonstrated to be suitable for RTHS. 412 5. Conclusions 413 In this study, a novel tracking controller for dynamics compensation in real-time 414 hybrid simulation was proposed. Although real-time hybrid simulation is an eﬀective 415 testing technique, it is challenging as the inherent dynamics of the actuation system 416 introducetimedelaysthatrequirecompensation. Therefore,agoodtrackingcontrolleris 417 anelementofthephysicalspecimentestingsetupneededtoobtainhighﬁdelityreal-time 418 hybridsimulation. Theproposedcontrolschemeisanextensionofthemodelpredictive 419 control presented in [21] to include an adaptive behavior. It is based on adaptive 420 model predictive control, a linear time-varying Kalman ﬁlter, and a real-time model 421 identiﬁcation algorithm. For the latter, ARX polynomial models are identiﬁed in real- 422 timetoestimatethechangingplantdynamicsandareusedtoupdatethepredictionmodel 423 of the model predictive controller. A parametric virtual real-time hybrid simulation 424 case study is used to validate the performance and robustness of the proposed control 425 scheme. In particular, the plant dynamics were varying during the simulation in the 426 presence of disturbance and measurement noise, additive to its output. Additionally, 427 parameters of the hybrid model were chosen and assigned with random values via 428 prescribed probability distributions. A total of 200 samples were generated and 200 429 hybrid simulation evaluations were performed to assess whether the proposed control 430 scheme is robust and preserves the desired performance. Results indicate that the 431 controller can guarantee small tracking errors under uncertainties that may be present 432 during the hybrid simulation. As a result, the eﬀectiveness of the proposed control 433 methodology was demonstrated proving that the adaptive model predictive control is 434 a suitable approach for real-time hybrid simulation. Future work aims to deploy the 435 proposed tracking controller in a test bench as well as to examine how multiple model 436 predictive controllers could coordinate in real-time, controlling multiple individual 437 physical substructuresof ahybridmodel. 438 25 DRAFT 6. Data Availability Statement 439 All data and code that support the ﬁndings of this study are available from the 440 corresponding author upon reasonable request. 441 Acknowledgments 442 ThisprojecthasreceivedfundingfromtheEuropeanUnion’sHorizon2020research 443 and innovation programme under the Marie Skłodowska-Curie grant agreement No. 444 764547. Thesoleresponsibilityofthispublicationlieswiththeauthor(s). TheEuropean 445 Union is not responsible for any use that may be made of the information contained 446 herein. 447