CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. Provisional Patent Application Serial No. 63/261,342, filed September 17, 2021, which is incorporated herein by reference for all purposes.

STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH

[0002] N/A

BACKGROUND

[0003] Many modern systems are formed from a variety of technologies, often including electrical systems, mechanical systems, computer systems, and any of variety of other systems or subsystems. For example, many industrial systems include mechanical infrastructure sub-systems, electrical sub-systems for operating the industrial system, and a computer sub -system for monitoring and/or controlling the operation of the electrical sub-systems or components.

[0004] In an effort to connect and coordinate operation of conglomerations of subsystems, modeling has been employed. Furthermore, more recently, model updating has emerged as a tool for system identification, parameter estimation, damage identification, response reconstruction, and virtual sensing. Model updating, the process of inferring models from data, is prone to the adverse effects of modeling error, which is caused by simplification and idealization assumptions in the mathematical models. The prediction error is usually modeled as an independent Gaussian white noise process in a Bayesian model updating framework.

[0005] In the application of model updating, the unknown model parameters and/or input loads are estimated by minimizing the mismatch between the measured and model-predicted responses. This mismatch, referred to as the prediction error, encapsulates the measurement noise and the effects of model uncertainties, which include model parameter uncertainties and modeling error. Model parameter uncertainties can be reduced through the model updating process given favorable identifiability conditions. Modeling error is caused by mathematical idealizations, approximations, and simplifications in the numerical model. If not accounted for properly, modeling error can cause estimation bias resulting in incorrect and/or inaccurate model updating outcomes.

[0006] In the classic non-adaptive recursive Bayesian model updating formulation, also referred to as non-adaptive Bayesian filtering, the prediction error is often assumed to be an independent Gaussian white noise process (i.e., a stationary, zeromean, uncorrelated Gaussian random process] for mathematical simplicity. This assumption results in a zero-mean Gaussian probability distribution function [PDF] for the prediction error with time-invariant covariance matrix. This limiting assumption can be violated in real-world conditions, perhaps most commonly due to the effects of modeling error. Incorrect characterization of the prediction error statistics will affect model updating performance and can result in biased estimation or divergence of updating parameters. It can also result in incorrect uncertainty quantification of model parameters, i.e., although the parameters may be estimated with reasonable accuracy, their uncertainty bounds may not be realistic.

[0007] To address this issue, Bayesian inference methods for joint model parameters, input loads, and noise identification have been proposed for structural and mechanical engineering applications. One Bayesian method for estimation of the diagonal entries of the prediction error covariance matrix has been proposed where the estimated covariance matrix is then fed into an extended Kalman filter (EKF) for joint stateparameter estimation. Other methods have used a dual adaptive filtering method to handle the effects of modeling error that consisted of a Kalman filter to estimate the diagonal entries of the prediction error covariance matrix based on a covariancematching technique, and an unscented Kalman filter (UKF) to estimate the unknown model parameters.

[0008] The above-mentioned studies assume that the prediction error is uncorrelated in space and, therefore, noise identification is limited to the estimation of the diagonal entries of the prediction error covariance matrix. Other methods have relied upon an adaptive Kalman filter using two types of covariance-matching methods, i.e., forgetting factor method and moving window method to estimate the full covariance matrix of the prediction error jointly with the unknown model parameters. Other methods have combined the Kalman filter method with a covariance-matching method to estimate the full covariance matrix of the prediction error jointly with the state vector and unknown model parameters. Yet other methods have proposed a Bayesian method for joint estimation of the mean vector and covariance matrix of the prediction error. In this approach, the mean vector of the prediction error was assumed to have a Gaussian distribution, while an inverse-Wishart distribution was considered for the covariance matrix of the prediction error. First, the distributions were updated based on Bayesian inference and then mean estimates of the updated distributions were used in the UKF algorithm for joint parameter and state estimation. Importantly, this work only considers the mean vector and covariance matrix of prediction error as time-invariant.

[0009] Although these methods alleviated some of the limiting assumptions for the prediction error, the zero-mean Gaussian assumption still remains in addition to modeling errors. Thus, there remains a need for systems and methods that are not undermined by the adverse effects of modeling error, and prediction errors from Gaussian white noise processes.

SUMMARY OF THE DISCLOSURE

[0010] The present disclosure overcomes the aforementioned drawbacks by providing systems and methods for creating and using digital models of real-world systems. In one example, the systems and methods provided herein maybe used to create and/or use a quantifiable model of a real-world system formed of a variety of subsystems, and manage or quantify error or bias of the model. In one non-limiting example, a learning network may be used to jointly estimate model parameters and the statistical characteristics of the prediction error that includes the effects of modeling error and measurement noise. The learning network may be an adaptive recursive Bayesian inference framework. The prediction error may be of the form of a non-stationary Gaussian process with unknown and time-variant mean vector and covariance matrix to be estimated. This provides the ability to account for the effects of time-varying model uncertainties in the model updating process. In a non-limiting example, a digital twin may be modeled for complex, real-world systems using a Bayesian physics-based digital twin based on an adaptive recursive Bayesian inference method.

[0011] In one configuration, a method is provided for managing a physical system with a computer system. The method includes generating, with the computer system, a physics-based model of at least one physical asset in the physical system. The method also includes collecting sensor data from sensors configured to record a parameter or a response of the at least one physical asset for training the physics-based model. The method also includes subjecting the physics-based model and the sensor data to a data assimilation process to generate or identify model parameters and/or input loads from the sensor data for the at least one physical asset and quantify a prediction error. The method also includes updating the physics-based model based upon the model parameters and the quantified prediction error to generate a physics-based digital twin and directing operation of the physical system using the physics-based digital twin and the input loads. In some configurations, the physics-based digital twin may be used for virtual sensing and response reconstruction, which is an approach to guide operation of the physical system using the physics-based digital twin.

[0012] In one configuration, a system is provided for modeling a physical asset. The system includes a computer system configured to: generate a physics-based model of the physical asset in a physical system and collect sensor data from sensors configured to record a parameter or a response of the physical asset for training the physics-based model. In some configurations, the sensors may be configured to record dynamic response quantities of the physical asset for training the physics-based model. The computer system may also be configured to subject the physics-based model and the sensor data to a data assimilation process to generate model parameters and/or input loads from the sensor data for the physical asset and quantify a prediction error. The computer system may also be configured to update the physics-based model based upon the model parameters and the quantified prediction error to generate a physics-based digital twin. In some configurations, the computer system may be further configured to direct operation of the physical system based upon the digital twin and the input loads. In some configurations, the input loads may be unmeasured input loads.

[0013] The foregoing and other aspects and advantages of the present disclosure will appear from the following description. In the description, reference is made to the accompanying drawings that form a part hereof, and in which there is shown by way of illustration a preferred embodiment. This embodiment does not necessarily represent the full scope of the invention, however, and reference is therefore made to the claims and herein for interpreting the scope of the invention. Like reference numerals will be used to refer to like parts from Figure to Figure in the following description. BRIEF DESCRIPTION OF THE DRAWINGS

[0014] FIG. 1 is a flowchart of non-limiting example steps for a Bayesian physicsbased model or digital twin based on an adaptive recursive Bayesian inference method.

[0015] FIG. 2A is a diagram of a non-limiting example 3-story 1-bay steel moment frame.

[0016] FIG. 2B is a graph of non-limiting example ground acceleration time history of an earthquake.

[0017] FIG. 3 depicts graphs of non-limiting example estimated model parameters for methods of FIGS. 2A and 2B.

[0018] FIG. 4 is a graph of non-limiting example estimated prediction error mean of FIGS. 2A and 2B.

[0019] FIG. 5 is a graph of non-limiting example covariance of the prediction error of FIGS. 2A and 2B.

[0020] FIG. 6 is a graph of estimated absolute acceleration responses at each floor compared with the true/nominal counterparts of FIGS. 2A and 2B.

[0021] FIG. 7 is a graph of the moment-curvature response and the stress-strain response of FIGS. 2A and 2B.

[0022] FIG. 8 depicts graphs of non-limiting example estimated model parameters for methods of FIGS. 2A and 2B.

[0023] FIG. 9 is a graph of non-limiting example estimated prediction error mean of FIGS. 2A and 2B.

[0024] FIG. 10 is a graph of non-limiting example covariance of the prediction error of FIGS. 2A and 2B.

[0025] FIG. 11 is a block diagram of an example of a computer system that can perform the methods described in the present disclosure.

[0026] FIG. 12 is a block diagram of an example digital twin formed from a nonlimiting example physics-based model, measurement data, and Bayesian inference approach.

DETAILED DESCRIPTION

[0027] Systems and methods are provided for creating and using physics-based digital twins of a physical system with a Bayesian data assimilation process. In one example, the systems and methods provided herein may be used to create and/or use a quantifiable model of a real-world system formed of a variety of sub-systems, estimate the unknown/unmeasured input loads on the system, and create and/or manage or quantify error or bias of the model. In one non-limiting example, an adaptive recursive Bayesian inference framework may be used to jointly estimate model parameters, input loads and the statistical characteristics of the prediction error that includes the effects of modeling error and measurement noise. The prediction error maybe of the form ofa non- stationaiy Gaussian process with unknown and time-variant mean vector and covariance matrix to be estimated. This provides for accounting for the effects of time-varying model uncertainties in the model updating process. In a non-limiting example, a digital twin may be modeled for complex, real-world systems using a Bayesian physics-based digital twin based on an adaptive recursive Bayesian inference method.

[0028] In a non-limiting example, a digital twin may be modeled for an offshore wind turbines (OWTs) system using a Bayesian physics-based digital twin based on an adaptive recursive Bayesian inference method. A digital twin or updated model may be used for virtual sensing to reconstruct the unmeasured response components to manage a system, such as an OWT, in order to reduce operation and maintenance costs and improve maintenance efforts.

[0029] The prediction capabilities of data-driven models are only as good as their training dataset, because they do not inherently capture the mechanics-based response behavior of the OWT. The training dataset will necessarily depend on the OWT type, design, and site conditions (e.g., soil properties, foundation installation, wind/wave loadings, etc.) that can vary significantly from one OWT to another. The limited amount of available data from OWT farms, and the variety of design and operational configurations result in the lack of reliable training dataset for supervised learning techniques, which in turn can result in unsatisfactory predictive capabilities of these data-driven models. Specialized data-driven models have been developed and calibrated based on experimental and real-world data to predict the remaining life of components (e.g., drivetrain bearings, consumable components, and the like). Although useful, these methods lack a holistic system-level approach; they are localized and component-specific and lose sight of the overall system dynamics. Moreover, data-driven models still require the installation of various sensors at different locations, which results in the shortcomings highlighted earlier. These methods lack the physics- and mechanics-based principals that can be used to characterize the complex response behavior of OWTs. Creating a digital replicate based on the mechanics of the complex OWT system may significantly broaden and expand the predictive maintenance capabilities.

[0030] In a non-limiting example, a framework is provided for remote monitoring, diagnosis, and prognosis of OWT systems and components to guide wind farm management and maintenance. The technology framework may be a Bayesian physicsbased digital twin based on an adaptive recursive Bayesian inference method. Physicsand mechanics-based models of the system may be used, and the model may be trained using sensory data.

[0031] The physics-based digital twin framework can be applied to an offshore wind turbine (OWT) system and its components, including foundation & soil, substructure & tower, drivetrain, blades, and to monitor and/or predict the lifetime of the involved components, and guide management and maintenance strategies. Components may include drivetrain bearings, consumable components, and the like. The methods in accordance with the present disclosure can be extended to fixed-based as well as floating OWTs. While the applications to different OWT components are different, they may all be based on a physics-based digital twin and may utilize the same mathematical backbone. Each application may use specific development, verification, and validation processes.

[0032] In a non-limiting example, the systems and methods in accordance with the present disclosure may be verified numerically using a 3 degree-of-freedom nonlinear mass-spring-damper system representing a 3-story 1-bay nonlinear steel moment frame excited by an earthquake. Comparison of the results with those obtained from a classical non-adaptive recursive Bayesian model updating method shows the efficacy of the systems and methods in accordance with the present disclosure in estimation of the prediction error statistics and model parameters.

[0033] The recursive Bayesian inference framework provides a mathematical basis to formulate model updating in the presence of uncertainties and noise. In this framework, the prediction error may be modeled as a random process characterized by a joint probability distribution function (PDF). The prediction error can be a correlated, non-stationary, non-white-noise, and non-Gaussian random process.

[0034] A recursive Bayesian inference formulation may be used to jointly estimate the unknown model parameters, unknown input loads, and the statistical characteristics of the prediction error. In non-limiting examples, the statistical characteristics include mean vector and covariance matrix of the prediction error, and the like. The prediction error may be modeled as a non-stationary Gaussian random process with time-variant mean vector and covariance matrix to be estimated iteratively and jointly with the unknown model parameters. In a non-limiting example, the estimation problem may be solved using a two-step marginal maximum a posteriori (MAP) estimation approach.

[0035] Referring to FIG. 1, a flowchart of non-limiting example steps for managing a complex physical system using a model is shown. In some configurations, a physicsbased digital twin may be integrated with measured data based on an adaptive recursive Bayesian inference method. A physics-based model of a physical system, asset or a component of a system may be generated at step 102. The framework may include remote monitoring, diagnosis, and prognosis of a physical system or component. In a nonlimiting example, the system is an OWT system and components to guide wind farm management and maintenance. Physical assets or components of an OWT include, but are not limited to blades, blade pitch system, gearboxes, protective coatings to prevent corrosion, hardware, bearings, tower, rotor, rotor shaft, brakes, grease, electrical components, control software, electrical generator, power cabling, wind sensors/anemometer or wind vane, yaw system, nacelle, and the like. The technology framework may be a Bayesian physics-based digital twinning based on an adaptive recursive Bayesian inference method. Physics- and mechanics-based models of the system may be used, and the model may be trained using sensory data.

[0036] Sensor data from the physical asset or system may be collected at step 104. Sensors may monitor wearable or consumable parts or components in order to track maintenance or other needs. Non-limiting examples include monitoring the wear of bearings, monitoring life expectancy of consumables (e.g. drivetrain components, and the like), determining the output of electrical generators, or assessing the response time of a yaw control system, the pitch angles a blade can achieve and the like. Sensors may include cameras, OPTICAL detectors, friction gauges, force gauges, accelerometers, strain gauges, clearance gauges, electrical current detectors, and the like. It is to be appreciated that any number of sensors may be used to monitor any number of systems, such as using a camera system along with a computer vision algorithm to measure the dynamic vibration or static deflection of a component, or to monitor the yaw control system. Any combination of sensors may also be used, such as to use a strain gauge to assess strain in a component in combination with an accelerometer to assess dynamic vibrations. The physics-based model may be trained with the sensor data at step 106.

[0037] Digital twin parameters for the physics-based model and errors of the model may be estimated at step 108, such as with a learning network or data assimilation process. In a non-limiting example, the learning network or data assimilation process uses a recursive Bayesian inference method. Digital twin parameters may include mechanical properties of the components, input dynamic forces, and the like, as described for sensors and monitoring above. The physics-based model may be updated based upon the estimated parameters and errors at step 110. A parameter of the physical system may be predicted, or operation of the physical system may be directed based upon the updated physics-based model representing a digital twin of the physical system at step 112. The predicted parameter or directed operation of the physical system may include virtual sensing, demand estimation, effects of wear, life prediction of the components, necessary maintenance, and the like.

[0038] In a non-limiting example, a physics- or mechanics-based model of the physical system may be utilized. Sensors may collect sparse data from the physical system. Through the framework, the information contained in data are integrated with the physics-based models using an adaptive Bayesian inference approach to reduce model uncertainties and provide improved prediction capabilities by estimating the modeling errors. The model, trained with data, represents a replicate of the physical system and is utilized for response and demand prediction, damage diagnosis, and remaining life prediction in various components.

[0039] In some configurations, the digital twin can account for the modeling errors and uncertainties, which in turn will be considered for estimation of remaining useful life of the OWT. The method uses the measured data to estimate jointly the input loads applied on the OWT and the parameters characterizing the physics-based model behavior. Through a dynamic sub -structuring approach, the complex OWT system is broken mathematically into different components, each may be characterized by a separate physics-based model. The physics-based model of each component will be trained with sensor data, to estimate the digital twin parameters and the error (both bias and variance) in a recursive manner as continuous monitoring data is collected. The methodology may underline the value of information and how having more sensors and data can reduce the estimation uncertainty of the remaining life prediction. [0040] Bayesian Inference Formulation for Joint Input, Model, and Noise Identification

[0041] A sequential window-based Bayesian filtering approach may be used for joint input, model and noise estimation using output-only vibration measurements such as acceleration, strain, velocity, or displacement. This approach may be applicable to linear and nonlinear structural or mechanical systems such as wind turbines, buildings, bridges, or aircrafts subjected to unknown (not measured) dynamic loadings such as wind, wave, earthquake, traffic, and the like. In this framework, the estimation time interval may be subdivided into N _w successive overlapping estimation. The estimation problem may be solved sequentially over sliding time windows.

[0042] The measured response vector of the structural system at estimation time window m, can be related to the response predicted through a numerical model, in which is an extended unknown parameter vector at the mth estimation window including the time-invariant unknown model parameter vector and the unknown input time history is ^e number of unknown model parameters, n _t is the number of unknown components in input time histories, and is the number of time steps in estimation window. is the vector of measured responses between time steps indicating the number of measurement channels. is the response function of the linear or nonlinear model (e.g., finite element) between time steps to an input force time history from time step . For the sake of notation brevity, is replaced by henceforth, denotes random vector x following a Gaussian (or Normal) distribution with mean vector x and covariance matrix Σ, and the PDF of X is shown as denotes prediction error vector accounting for the mismatch between the measured and model-predicted responses of the structure and may be modeled as a non-stationary Gaussian random process with unknown and time-variant mean vector and covariance matrix to be estimated recursively and jointly with the unknown model parameters and input time history .

[0043] The prediction error may be independent across time. Furthermore, the mean vector and covariance matrix of the prediction error may be time invariant at each estimation window, i.e., for all where are time-invariant mean vector and covariance matrix of the prediction error at mth estimation window, respectively. Therefore, the estimation problem at each estimation window may be limited to estimation of and .

[0044] To proceed with the Bayesian inference formulation may be modeled as random vectors and may be modeled as a random matrix. The estimation problem may be solved sequentially over all the estimation windows, and at each estimation window m, the framework may include "prediction” and "correction” steps. In the prediction step, the joint distribution of updating parameters, are propagated from the previous window to the current estimation window through a dynamic model and they are used as prior information. In the correction step, prior estimates of the updating parameters, denoted as are corrected to posterior estimates, denoted as using an iterative process through the Bayes’ theorem by absorbing the information in the measurement between time steps and

[0045] Bayesian Inference Formulation

[0046] Using the Bayes’ theorem, the joint posterior distribution of updating parameters, given the measured responses from time steps , can be derived as where is the likelihood function and is the joint prior distribution of at mth estimation window. The normalizing (evidence] term may be ignored in Eq. (2); therefore, the sign denoting "proportional to” is used. The terms in the right-hand side of this equation may be expanded. Based on Eq. (1), the likelihood function follows a Gaussian distribution as where based on the assumption that the perdition error is independent across time.

[0047] The joint prior distribution, , can be written in a hierarchical form as in which is referred to as hyper-prior with hyper-parameters of By substituting Eq. (4) into Eq. (2), it can be followed that

[0048] The prior distribution of may be approximated as Gaussian, i.e.,

[0049] where mean vector is the prior estimate for is the prior covariance matrix of . Furthermore, the prior distribution of may follow a Normal-Inverse-Wishart [NIW) distribution. In Bayesian statistics, NIW distribution may be used as the joint conjugate prior for the mean vector and covariance matrix of a Gaussian distribution. The conjugacy guarantees the same functional form for the posterior and prior distributions. NIW distribution is the product of a Normal [or Gaussian) distribution and an Inverse-Wishart (IW) distribution. The prior distribution of in Eq. [5] can be expressed as where is the prior mean vector of [also considered as prior estimate for is the confidence parameter, is the degree of freedom parameter and is a symmetric positive definite scale matrix. Considering that , the right-hand side of Eq. [7) can be separated as follows. [0050] The Normal and IW distributions in Eq. (8-a) and (8-b) are defined below, ignoring normalizing terms.

The terms | . | and tr(.) in these equations denote matrix determinant and matrix trace, respectively.

[0051] Dynamic Models

[0052] The evolution of updating parameters from estimation window m — 1 to m maybe characterized by a series of dynamic models. For propagating model parameter in the extended parameter vector , we can have For the input load time history, the estimated unknown input time history over the estimation window may be divided into two parts . The first part, from , which has overlap with the previous estimation window, , can be transferred from the previous window and used as the initial estimate for the current estimation sequence. The second part, from can be initialized with the average value of the first part. Therefore, can be written as where t _s is the sliding (or moving) rate defined as the difference (in number of time steps) between the starting point of two consecutive windows, i.e., ₁ is a vector of ones with the length of t _s .

[0053] The uncertainties associated with model parameters and the first part of input load time history can be transferred from the previous window to the current window by deriving the conditional posterior covariance matrix from the previous window. The covariance matrix for the second part of the estimated input load time history can initialized as a diagonal matrix with identical diagonal entries. Therefore, we have where the conditional posterior covariance matrix maybe defined as

[0054] To improve the convergence of the iterative procedure, a disturbance Q may be added to the posterior covariance matrix at each iteration to provide the prior covariance matrix. is a constant diagonal matrix with small positive diagonal entries.

[0055] An explicit presentation of a dynamic model for the covariance matrix of the prediction error, may not be available. The following two dynamic models for the statistical parameters of the prior IW distribution may be used, based on a heuristic model. where p ∈ ((,1] is a forgetting factor. If p = 1, Eqs. (13-a) and (13-b) results in a stationary prediction error covariance matrix, while smaller values of p allows for larger time variations in the statistical properties of covariance matrix. The mode (most probable) value of the prior IW distribution, considered as prior estimate for , is defined as while the posterior estimate for is computed as Based on the dynamic models defined in Eq. (13), it can be concluded that .

[0056] Similar to the dynamic model considered for a heuristic model may be considered for propagating the uncertainties of the through time as where is another forgetting factor. Similar to the case presented for represents a stationary , while smaller values of p' result in larger time variations in the statistical properties of prediction error mean.

[0057] Two-step Marginal MAP Estimation

[0058] Posterior estimates of updating parameters may be found using a maximum a-posteriori [MAP] estimation method. MAP corresponds to the parameters that maximize the joint posterior distribution

[0059] The joint posterior distribution can be factorized into two marginal distributions as follows using the product rule.

[0060] Based on the marginal MAP estimation method, the two terms on the right- hand side of Eq. (16) may be maximized separately to find the MAP estimates, i.e.,

[0061] The correction step in the Bayesian inference formulation may be divided into two separate MAP estimation problems that can be solved iteratively to converge to the MAP estimates of the joint posterior distribution , as described below.

[0062] Marginal MAP Estimate of

[0063] The MAP estimation problem in Eq. (17) may be used in non-adaptive Bayesian model updating formulations. Based on the Bayes’ theorem, it can be written that

[0064] With both the likelihood and the prior distributions for being Gaussian according to Eqs. (3) and (6), respectively, the posterior distribution will also be Gaussian, i.e., . The mean vector (considered the same as the MAP estimate for ), and the covariance matrix can be obtained as:

The term is the model response sensitivity matrix with respect to . The terms are the posterior estimate of the mean vector and covariance matrix of prediction error. Based on the assumption that they are time-invariant at each estimation window, they can obtain as where are the MAP estimates of the time-invariant mean vector and covariance matrix of the prediction error at estimation window m; they may be estimated by solving Eq. (18) as described below. Eqs. (20)-(21) are similar to the equations used in the non-adaptive Bayesian model updating methods with the exception of which is now present in Eq. (20) to account for the non-zero-mean prediction error.

[0065] Marginal MAP Estimates of

[0066] From Eq. (18), the term can be written using the Bayes’ rule as

[0067] According to Eq. (7), the prior distribution has a NIW distribution, which is a conjugate prior for the Gaussian likelihood. Therefore, the posterior will also be with updated parameters

[0068] Based on how the perdition error may be independent across time and the mean vector and covariance matrix of the prediction error are time-invariant at each estimation window, by substituting Eqs. (3) and (7) into Eq. (23), the following updating equations can be derived. where is the average value of prediction error in the estimation window. The mode of the posterior NIW distribution may be selected as the MAP estimates for i.e., as shown in Eqs. (24-a) and (24-e).

[0069] The solution of the MAP estimation problem in Eq. (15) result in solving the coupled Eqs. (20), (24-a), and (24-e) simultaneously using a fixed-point iteration algorithm. First , may be estimated using Eq. (20) given the obtained from the previous iteration; then may be updated using Eqs. (24-a) and (24-e) based on the estimated , and these estimates may be iterated until convergence is achieved. In a non-limiting example, three convergence criteria can be considered based on relative L2 norm of the difference between two consecutive estimations of . A convergence tolerance may be used, such as a convergence tolerance of 0.01. This value can be adjusted by balancing accuracy versus the computational cost. An L2 norm of a matrix is equal to its largest singular value. Table 1 provides a non-limiting example flowchart of a two-step marginal MAP estimation approach.

[0070] Table 1. Two-step marginal MAP estimation approach for joint Input, model and noise identification.

[0071] Non-limiting example Verification: 3 -story 1-bay Steel Moment Frame

[0072] Model Description and Data Simulation

[0073] Referring to FIGS. 2A and 2B, a non-limiting example 3-story 1-bay steel moment frame is shown in FIG. 2A, and a ground acceleration time history of 1989 Loma Prieta earthquake recorded at Los Gatos station in the 0° component is shown in FIG. 2B. Performance of the Bayesian inference framework for joint estimation of model parameters and noise (noise or prediction error = modeling error + measurement noise) may be evaluated with the nonlinear model of a 3-story 1-bay steel moment frame under earthquake excitation as shown in FIGS. 2A and 2B. A 2D nonlinear FE model of the structure was developed in the open-source FE analysis platform OpenSees. Columns are made of A992 steel with W14x311cross-section and beams were made of A36 steel with W24x68 cross-section. Nodal mass of 80 metric tons was considered at beam-column nodes shown. Columns and beams were modeled using force-based beam-column elements with fiber sections. Seven integration points were considered for numerical integration along the length of each element using Gauss-Lobatto quadrature rule. Column and beam webs were discretized into 10 fibers along height and one fiber across width. Their flanges were also discretized into one fiber along height and 3 fibers along width. The uniaxial material model for steel fibers was based on the Giuffre-Menegotto- Pinto (GMP) constitutive model. Rayleigh damping was considered to model the structural damping assuming 2% damping ratio for the first two modes.

[0074] The 1989 Loma Prieta earthquake (0° component at Los Gatos station) was selected as the base excitation as illustrated. The horizontal absolute acceleration response time histories at each floor (referred to as true/nominal responses and denoted by y true ) were simulated and contaminated with measurement noise to represent the sensor measurements (denoted by y). The measurement locations are shown. The measurement noise was considered as a non-stationary Gaussian random process with time-variant mean vector and covariance matrix denoted by , respectively. The measurement noises were assumed statistically uncorrelated; therefore, the off- diagonal entries of at each time step k were equal to zero. Sine functions were considered to model variation of and the diagonal entries of in time. where N is the number of time steps, and are constant coefficient vectors and defined as follows.

(27)

[0075] Three parameters of the GMP model for beams and columns were considered as the unknown (updating) model parameters, including the Young’s modulus E, yield stress σy , and strain-hardening ratio b. The parameters were specified by subscripts b and c for beams and columns, respectively. The nominal (or true) values of these parameters used for simulation are reported in Table 2. The mass and stiffness proportional components of Rayleigh damping, i.e., a and β , respectively, were also considered as the unknown model parameters to be estimated. Their true or nominal values were considered as a ^true = 0.1718 s- ¹ and β ^true = 0.0014 s . Therefore, the unknown model parameter vector 0 included 8 parameters, each one was normalized by its true value, i.e.,

[0076] Table 2. True/nominal values of the three parameters of GMP model for columns and beams used for measurement simulations.

[0077] Identification Results

[0078] The methods in accordance with the present disclosure were applied to estimate the unknown model parameter vector θ together with the mean vector and covariance matrix of prediction error. The initial value of θ was selected as = [0.7 0.7

1.2 0.8 1.3 0.8 1.2 0.7] ^T with the initial covariance matrix , i.e., a 20% initial coefficient of variation (COV) was considered to characterize the prior covariance

, , , matrix. The process noise covariance matrix was assumed as nd the forgetting factor parameters are selected as p = 0.95 and p' = 0.95. The initial mean vector and covariance matrix of the prediction error were assumed as and respectively, where vector is defined in Eq. (26). Other initial parameters considered for the NIW distribution are , and . In this verification study, the results of the proposed joint model and noise identification method is compared with a classic non-adaptive recursive Bayesian method (Ebrahimian et al. 2015), in which the prediction error is assumed to be an independent zero-mean Gaussian white noise with a time-invariant , diagonal covariance matrix, i.e.,

[0079] Referring to FIG. 3, graphs of non-limiting example estimated model parameters for both methods of FIGS. 2A and 2B are shown. The Bayesian method can estimate accurately all 8 unknown model parameters, while model parameter estimates are biased or incorrect for the non-adaptive Bayesian model updating method. The non- adaptive Bayesian model updating method can even result in divergence of the model updating process or producing nonphysical model parameter estimates, e.g., zero estimated values for Eb.

[0080] Referring to FIG. 4 and FIG. 5, graphs of non-limiting example estimated mean and covariance of the prediction error are shown, respectively, at each time step. It can be seen that the Bayesian method accurately tracks the trend of the true mean vector and covariance matrix in time.

[0081] Referring to FIG. 6, graphs of estimated absolute acceleration responses at each floor - using the final estimated model parameters - are compared with the true/nominal counterparts. The noticeable discrepancies between the estimated and nominal responses for the non-adaptive Bayesian model updating method, which is due to the biased model parameter estimates, clearly show the incapability of the non- adaptive Bayesian model updating method to perform in the presence of a time-varying prediction error. The acceleration responses predicted using the joint model and noise identification method match the nominal responses well.

[0082] Referring to FIG. 7, graphs of the moment-curvature response at the base section of the left column - Section 1-1 in FIG. 2A - and the stress-strain response of the top flange of the first-floor beam at the plastic hinge location - Section 2-2 in FIG. 2A - are shown and were estimated from the updated models using the non-adaptive and adaptive model updating methods. The detrimental effects of the time-varying noise in the response prediction and virtual sensing capability of the non-adaptive Bayesian model updating may be seen. The adaptive Bayesian model can handle the time-varying noise effects and provide accurate virtual sensing capability.

[0083] Non-limiting examples have been provided of an adaptive recursive Bayesian inference framework for joint model and noise identification using a two-step marginal maximum a posteriori [MAP] estimation approach. Noise or prediction error can include the measurement noise and the effects of modeling error. Non-limiting example separate MAP estimation problems may be solved iteratively: one to estimate the unknown model parameters and another to estimate the mean vector and covariance matrix of the prediction error. Verification of a non-limiting example included using numerically simulated data obtained from a nonlinear steel moment frame model subjected to earthquake excitation. The absolute acceleration responses at each floor were simulated and polluted by a non-stationaiy Gaussian noise with time-variant mean vector and covariance matrix. Eight model parameters, including six parameters characterizing the constitutive models of the beams and columns steel material, and two Rayleigh damping coefficients, were considered as unknown and estimated using the proposed approach and a non-adaptive recursive Bayesian model updating approach for comparison. The non-limiting example verification demonstrated the detrimental effects of time-varying noise on the non-adaptive Bayesian model updating results, where the estimated model parameters were significantly biased.

[0084] An adaptive or recursive Bayesian inference framework for joint model and noise identification was able to estimate the model parameters correctly, due to its capability to estimate the time-variant mean and covariance of the prediction error. Considering the statistical characteristics of prediction error as unknowns to be estimated provides additional degrees of freedom in the recursive Bayesian model updating approach and can alleviate the biased or incorrect model parameter estimation results. The method can provide a step-forward to account for the effects of modeling error in finite element model updating.

[0085] Another non-limiting example can include estimating both the unknown model parameters and statistical characteristics (mean vector and covariance matrix] of the prediction error, as well an approximation of their joint posterior distribution. Exact calculation of this high-dimensional joint posterior distribution is intractable, so the process requires approximation. Two approximation schemes can be used: stochastic or sampling methods such as, for example, Markov chain Monte Carlo (MCMCJ, and deterministic or variational frameworks such as, for example, variational Bayesian (VB). In comparison to the sampling methods, the VB method is analytically tractable and is computationally less demanding. The VB method has been successfully applied for joint state and noise estimation in navigation, target tracking, and control related applications. In these applications, adaptive VB Kalman filter method was used to estimate jointly the covariance matrix of a zero-mean prediction error and the state of linear or nonlinear state-space models. The VB method assumes that the approximate joint distribution is the product of some single- or multi- variable factors and uses the Kullback-Leibler (KL) divergence to minimize the difference between the approximation and the true posterior. [0086] Continuing, another non-limiting example can provide an adaptive method for nonlinear model updating based on the VB method to approximate the joint posterior distribution of the unknown model parameters and statistical characteristics of the prediction error at each time step.

[0087] Below, a model updating problem statement and a detailed derivation of the proposed VB method for estimating the joint posterior distribution of unknown model parameters and statistical characteristics of the prediction error are provided. Further, the formulation of the VB method is then compared with the two-step marginal MAP estimation method, and then the proposed method is verified using a 3- stoiy 1-bay nonlinear steel moment frame under earthquake excitation, and the results are compared to those from the two-step marginal MAP estimation method provided above and a non-adaptive Bayesian model updating method.

[0088] Model Updating Problem Statement

[0089] Consider the measured response of a nonlinear (or linear] dynamic system y, and its corresponding model (e.g., finite element (FE)) prediction h(0) where 0 is the vector of unknown model parameters. The parameter estimation problem at time k = 1, 2, ... , N can be formulated as where is the process noise and is the prediction error. The input forces are assumed to be known, so for the sake of notation brevity the dependency of model prediction response to the input forces is not shown explicitly in Eq. (2). The process noise is assumed to follow a zero-mean Gaussian white noise process with covariance matrix In non-adaptive Bayesian model updating methods, the prediction error is assumed as a zero-mean Gaussian white noise process with a constant or time-invariant covariance matrix R, i.e., For the parameter estimation problem defined in Eqs. [1] and (2), the non-adaptive Bayesian methods can be used to find an estimate for the first two statistical moments of unknown model parameters. However, in the adaptive Bayesian model updating methods, the prediction error can be modeled as a non-stationary Gaussian random process with unknown and time-variant mean vector and covariance matrix , to be estimated recursively and jointly with the unknown vector of model parameters As discussed above, a two- step maximum a posteriori (MAP) estimation method was used to estimate by maximizing the joint posterior distribution , i.e.,

[0090] To solve this MAP problem, the problem was broken into two iterative MAP estimation problems as

An example aim of this method is to find the whole joint posterior distribution of unknown model parameters and noise . Nevertheless, analytical working with this joint posterior distribution is difficult because the number of variables is high, and the joint distribution is highly complex. To overcome this issue, the VB method is used to approximate the joint posterior distribution. The VB method and the derivation details are provided below.

[0091] Variational Bayesian (VB) Method

[0092] VB is a method to approximate a joint distribution p by a joint distribution Q which can be factorized into single-variable or grouped-variables factors. The Kullback-Leibler (KL) divergence criterion is then used to make Q as close as possible to p. Using the VB method, the joint posterior distribution of the unknown model parameters and prediction error mean vector and covariance matrix are estimated by separating the joint posterior distribution into two factors as follows, where are unknown distributions which can be obtained by minimizing the KL divergence between the right- and left-hand side of Eq. (4). The KL divergence is defined as

[0093] Using variational calculus to minimize the above KL-divergence with respect to each of while keeping the other one fixed will result in Eqs. (6) and (7). where denotes the expected value of f(x) with respect to x with probability density function of . The terms denote the constants with respect to variables , respectively. Since Eqs. (6) and (7) are coupled through the term , analytical solutions are not available. Therefore, the fixed-point iteration algorithm can be employed to find approximate solutions for Eqs. [6] and [7]. To this end, the right-hand sides of Eqs. (6) and (7) are expanded from their innermost parentheses to the outer through the following steps.

[0094] The joint distribution in Eqs. (6) and (7) can be factored as, where is the likelihood function, is the prior distribution of is the prior joint distribution of is known because it depends on the past measurements.

[0095] Based on Eq. (2), and further expansion of the terms on the right-hand-side of Eq. (8), the likelihood function has a Gaussian distribution as

[0096] For the second and third terms on the right-hand-side of Eq. (8), it is assumed, similar to method described above that have prior distributions of Gaussian and Normal-Inverse-Wishart (NIW), respectively. NIW distribution is the product of a Gaussian (or Normal) distribution and an Inverse- Wishart (IW). The NIW is selected for the prior distribution because it is a conjugate prior for a Gaussian likelihood with unknown mean vector and covariance matrix. The conjugacy guarantees the same functional form for the posterior and prior distributions (O'Hagan and Forster 2004). Therefore, the second and third terms on the right-hand-side of Eq. (8] can be written as follows. where are the mean vector and covariance matrix of the unknown model parameters given measurements but not The minus superscripts represent the prior estimates are the prior estimates of statistical parameters used in the NIW distribution is the prior estimate for is symmetric positive definite scale matrix. are the confidence parameter and the degree of freedom parameter, respectively. They are scalar parameters and shall satisfy 1 where n _y is the number of measurement sensors.

[0097] By substituting Eqs. [9], (10) and (11) into Eq. (8), it can be followed that

[0098] The Gaussian (or Normal) distribution and the Inverse- Wishart (IW) distribution are proportional to the following expressions. where |. | represents determinant and tr(. ) denotes trace of a matrix. The sign representing "proportional to” is used in Eq. (13), as the normalizing terms are ignored. [0099] Using Eq. (12) and the definitions of Normal and IW distributions in Eq. (13), the term , which is used in Eqs. (6) and (7), can be expanded as follows. where denotes a constant with respect to all variables Now, having the expansion of ln in Eq. [14], the expectation terms in Eqs. [6] and [7] can be calculated. By getting the expectation from each term of Eq. (14), the expectation term in Eq. (6), can be expressed as follows. where ce is the summation of expectations of all terms in the right-hand side of Eq. (14) except the second and fourth terms and is constant with respect to . Note that the expectation of the fourth term of Eq. (14) is equal to itself because it does not depend on pt, and Rfc.

[00100] Using Lemma 1 in the Appendix (provided below), the expectation term in Eq. (15) can be evaluated as

[00101] Now, consider as a NIW distribution, i.e., where the plus superscrips represent the posterior estimates. Therefore, by substituting the mean and covariance of ), and taking the expectation with respect to Rk, the equation becomes [00102] It should be noted that in deriving Eq. which is the mean of IW distribution denoted as and can be calculated as follows.

[00103] Looking back to Eq. (6), is proportional to the exponential of So, by substituting Eq. (17) into Eq. (15) and taking exponential of Eq. (15), it can be followed that

[00104] By linearizing in Eq. (19) by the first-order Taylor expansion about is the sensitivity matrix of FE model with respect to

[00105] The first exponential term in the right-hand side of Eq. (20) shows a Gaussian distribution for ignoring the normalization term. By using Lemma 2 in the Appendix, the second exponential term in the right-hand side of Eq. (20) also represents a Gaussian distribution for 0t. Therefore, the right-hand side of Eq. (20) is the product of two Gaussian distributions which based on Lemma 3 in the Appendix, results in a Gaussian distribution, i.e., To obtain the parameters of this Gaussian distribution, we match the terms in the left-hand and the right-hand side of in Eq. (20), which results in the following two equations.

[00106] In a similar way, we can evaluate the expectation term in Eq. (7) as follows. Getting the mathematical expectation of Eq. (14) with respect to leadsto

where is sum of the expectation of the third, fourth, and the last term of the right-hand side of Eq. (14), and is constant with respect to

[00107] Now, consider and linearize by the first- order Taylor expansion about is the sensitivity matrix of the model with respect to , the expectation term in the right- hand side of Eq. (22) can be obtained as follows using Lemma 1 in Appendix.

[001 08] Based on Eq. (7), is proportional to the exponential of ). Substituting Eq. (23) into Eq. (22) and taking exponential of Eq. (22), can be found as follows.

[00109] The right-hand side of Eq. (24) includes the product of two Gaussian distributions for and one IW distribution for . The product of these two Gaussian distributions leads to a scaled Gaussian distribution based on Lemma 3 in the Appendix. Therefore, the right-hand side of Eq. (24) results in a NIW distribution, which is a product of a Normal distribution for and an IW distribution for By substituting in the left-hand side of Eq. [24], and using Lemma 4 in the

Appendix, the four parameters of NIW distribution can be derived as follows by matching the terms in the left- and right-hand sides of Eq. (24).

[00110] Now having as an approximation for the expectation of can be evaluated to represent point estimates of unknown model parameters and noise. Therefore, represented in Eqs. (21-a), (25-a), and (18) can be considered as the point estimates for Eqs. (21-a), (25-a), and (18) are coupled equations that can be solved iteratively using a fixed-point iteration algorithm. It is worth noting that Eqs. (21-a) and (21-b) are similar to non- adaptive Bayesian model updating formulations except the term which is added in Eq. (21-a) to consider non-zeromean prediction error.

[00111] The proposed algorithm of recursive VB for joint model and noise identification is presented in Table 3 below. Like other recursive Bayesian model updating algorithms, the proposed algorithm has two steps at each time : “prediction" and "correction." In the “prediction" step, the new measurement at time k is not given to the estimation process yet. Therefore, the prior estimates of the statistical parameters of the Gaussian distribution of , and the NIW distribution of , denoted by minus superscript, are predicted through a dynamic model using their posterior estimates at the previous time step k-1. Eq. (1) can be used as the dynamic model for unknown model parameters of to obtain . For predicting the prior estimates of the four parameters of NIW distribution, the dynamic models defined above can be used considering the forgetting factor parameters of . In the "correction" step, the prior estimates are corrected by the new measurement yt to get the posterior estimates, denoted by plus superscript, using Eqs. (21), (25), and (18) in an iterative way. In each iteration, the model-predicted responses and the sensitivity matrix need to be updated based on the updated unknown model parameters 0^ . However, calculating sensitivity matrix at each time step can be computationally demanding. To reduce the execution time, the prior sensitivity matrix in Eq. (25-d) can be used. Since the convergence criteria are not changed, using instead of has no effect on the final estimation results.

[00112] Table 3. Algorithm for recursive variational Bayesian (VB) method. [00113] Comparison of the VB Method with Two-step Marginal MAP Estimation Method

[00114] In this section, the proposed Variational Bayesian (VB) method is compared with the two-step marginal MAP estimation method discussed above for joint estimation of unknown model parameters and the mean vector and covariance matrix of the prediction error. The formulations of the two methods are closely similar. There are two main differences between the two methods as explained below. First, in the two- step marginal MAP estimation method the mode of IW distribution is used as point estimate, while in the VB method, the mean of IW distribution is assigned as (as shown in Table 4 below). The reason for this difference comes from the fact that the MAP approach is used in the former method, while the expectation value is used in the latter. It is worth noting that the mode and mean of IW distribution are not coincident. Second, different equations are used to calculate the term as shown in Table 4. The equation used to calculate in the VB method has an additional term, i.e.,

[00115] In most applications, these two differences have negligible effects on the results because of the following reasons. First, the difference between the mode and mean of IW distribution decreases in higher time steps as the value of increases in time in the denominator of the characterizing equations of Second, as the covariance matrix of the unknown model parameters decrease through the Bayesian estimation, the effects of this extra term on the results can also become negligible. Therefore, in most applications the final estimation results from the two methods are similar.

[00116] Table 4. Differences between the proposed VB method and two-step marginal MAP estimation method.

[00117] Numerical Study

[00118] In this section, the estimation results of the two disclosed methods are compared through a numerical study. Performance of the proposed VB method is similarly evaluated when applied to a numerical model of a 3-story 1-bay steel moment frame structure under earthquake excitation. The estimation results are compared with the two-step marginal MAP estimation method discussed above and a non-adaptive Bayesian model updating method. The story height and the bay width are 3.5 m and 6.0 m, respectively, as shown in FIG. 2A. The frame geometric and material properties are similar to the one used for the two-step marginal MAP estimation method. The numerical model of the frame is developed in OpenSees. For this, force-based beam-column elements with seven integration points are used for columns and beams. A single fiber is used to represent each flange of beam and column cross sections, while ten fibers are used to discretize their webs. The uniaxial Giuffre-Menegotto-Pinto (GMP) material model with primary parameters is used to model the steel fibers and simulate the nominal/true dynamic response of the structure, where E = Young’s modulus, Fy = yield stress, and b = strain-hardening ratio. The first three parameters denoted by subscript "c" are for columns, and the last three parameters denoted by subscript “b” are used for beams. A nodal mass m = 80,000 kg , shown by black circle in FIG. 2A, is considered for each story at the beam-column nodes to represent dead and live mass. T o model damping energy dissipation, Rayleigh damping with 2% damping ratio is considered for the first two vibration modes of the structure. To simulate measurement data, the frame structure is excited by Loma Prieta earthquake (0° component at Los Gatos station ] as shown in FIG. 2B. Then, the horizontal absolute acceleration response time histories of each floor (shown by black boxes in FIG. 2A] are extracted and contaminated with artificial measurement noise to result in simulated measurement data. The measurement noise is considered as a non-stationary Gaussian random process with time-variant mean vector and covariance matrix as follows.

[00119] The goal is to estimate the unknown model parameters and compare them with their true values The initial estimate of the unknown model parameters and its covariance matrix are selected as respectively. The initial mean vector and covariance matrix of the prediction error are assumed as respectively. Other initial parameters of the NIW distribution are selected as _y with n _y - 3. The process noise covariance matrix is selected as and the forgetting factor parameters used for defining dynamic model of the mean vector and covariance matrix of the prediction error are assumed as p = 0.95 and p' = 0.95. Now, the proposed VB method is applied to estimate jointly the unknown model parameter vector 0 and the mean vector and covariance matrix of prediction error. In this verification study, the results are compared with the two-step marginal MAP estimation method discussed above and the non-adaptive Bayesian model updating method when using the same initial values. As mentioned before, in the non-adaptive Bayesian method, a zeromean Gaussian white noise with a time-variant diagonal covariance matrix is assumed for the prediction error

[00120] This non-limiting method example use the VB approach and proposes an adaptive VB model updating method for joint model and noise identification. Detailed mathematical derivation is provided herein. To evaluate the performance of the proposed method, a 3-story 1-bay nonlinear steel moment frame subjected to earthquake excitation was used, in which six parameters characterizing the constitutive models of the steel beams and columns were considered as unknown. Absolute acceleration responses at each floor contaminated by Gaussian noise with time-variant mean vector and covariance matrix were considered as measurement data. The estimation results showed that the proposed VB method can accurately estimate both unknown model parameters and time-variant mean vector and covariance matrix of the prediction error. The proposed VB method was also compared to a two-step marginal MAP estimation method, and the non-adaptive Bayesian model updating method. The results showed that both adaptive methods have comparable performance while the non-adaptive method resulted in significantly biased estimations due to the adverse effects of non-stationary measurement noise.

[00121] FIG. 8 shows the time histories of the unknown model parameters estimated by all three methods: the proposed VB method, the two-step marginal MAP estimation method, and the non-adaptive Bayesian method. It can be observed that both the proposed VB method and the two-step marginal MAP estimation method can accurately estimate all six unknown model parameters. However, the non-adaptive model updating method converges to incorrect unknown model parameters, or even diverges. This shows that estimation of the statistical parameters of prediction error affects the model updating results significantly when modeling error is time-variant and therefore, shall not be ignored.

[00122] The estimated mean and variance of the prediction error (measurement noise in this example) for each channel are compared with their true values in FIG. 9 and FIG. 10, respectively. The non-adaptive Bayesian model updating method, as mentioned before, cannot estimate the mean and variance of the prediction error, so they remain constant during the estimation process. However, both the VB and the two-step marginal MAP estimation methods can accurately track the trend of the true/nominal mean and covariance of error through the time. There is a negligible difference between the estimations obtained from these two adaptive methods.

[00123] Referring now to FIG. 11, a block diagram of an example of a computer system 800 that can perform the methods described in the present disclosure is shown. The computer system 800 generally includes an input 802, at least one hardware processor 804, a memory 806, and an output 808. Thus, the computer system 800 is generally implemented with a hardware processor 804 and a memory 806. A real-world system 801, such as a physical system, may include a component or physical asset 803 as a subsystem.

[00124] In some embodiments, the computer system 800 can be a remote server or cloud-based system. The computer system 800 may also be implemented, in some examples, by a workstation, a notebook computer, a tablet device, a mobile device, a multimedia device, a network server, a mainframe, one or more controllers, one or more microcontrollers, or any other general-purpose or application-specific computing device. [00125] The computer system 800 may operate autonomously or semi- autonomously, or may read executable software instructions from the memory 806 or a computer-readable medium (e.g., a hard drive, a CD-ROM, flash memory), or may receive instructions via the input 802 from a user, or any another source logically connected to a computer or device, such as another networked computer or server. Thus, in some embodiments, the computer system 800 can also include any suitable device for reading computer-readable storage media. [00126] In general, the computer system 800 is programmed or otherwise configured to implement the methods and algorithms described in the present disclosure. For instance, the computer system 800 can be programmed to implement a learning network or data assimilation process, such as an adaptive or recursive Bayesian inference framework for digital twin modeling.

[00127] The input 802 may take any suitable shape or form, as desired, for operation of the computer system 800, including the ability for selecting, entering, or otherwise specifying parameters consistent with performing tasks, processing data, or operating the computer system 800. In some aspects, the input 802 may be configured to receive data, such as sensor data from the real-world system 801 where sensors may be used to monitor the real-world system 801 or the physical asset 803 as described above. In a non-limiting example, the real-world system is an OWT, and the physical asset is a component of the OWT. Such data may be processed as described above to train a physics-based model. In addition, the input 802 may also be configured to receive any other data or information considered useful for an adaptive or recursive Bayesian inference framework using the methods described above.

[00128] Among the processing tasks for operating the computer system 800, the one or more hardware processors 804 may also be configured to carry out any number of post-processing steps on data received by way of the input 802.

[00129] The memory 806 may contain software 810 and data 812, such as data acquired with the sensors, and may be configured for storage and retrieval of processed information, instructions, and data to be processed by the one or more hardware processors 804. In some aspects, the software 810 may contain instructions directed to a program to implement an adaptive or recursive Bayesian inference framework.

[00130] In addition, the output 808 may take any shape or form, as desired, and maybe configured for displaying results or a prediction for a physical system based upon the results of the physics-based model in addition to other desired information.

[00131] Referring to FIG. 12, a block diagram is shown for generating a physics based digital twin 908 from physics based model 902, measurement data 904, and Bayesian inference approach 906. The physics based model 902 characterizes the dynamics of the system as described above, such nonlinear behavior of foundation, dynamic interaction of the soil and support structure, and the like. In some configurations, the physics based model 902 may be a finite element model based on physically-meaningful parameters, such as foundation stiffness, aerodynamic damping parameters, and the like, as non-limiting examples. The exact values of the finite element models may be unknown or uncertain input exciting forces, such as wind and wave loads that are unknown and unmeasurable. Inherent modeling errors that stem from the mathematical idealizations and approximations may be included in the physics base model 902.

[00132] In a non-limiting example, measurement data 904 may include heterogeneous measurements on an OWT using sparse set of sensors. Bayesian inference approach 906 may include an estimate model input parameters, and or may characterize modeling errors.

[00133] The physics based digital twin 908 may be calibrated using estimated input parameters. As described above, the physics based digital twin 908 may estimate the consumed fatigue life in different components.

[00134] Appendix

[00135] This appendix includes four Lemmas used in deriving VB method relationships as discussed above.

[00136] Lemma 1: Expectation of a quadratic form

[00137] Considering vector , and symmetric matrix , the expectation of quadratic form can be obtained as follows where tr(. ) denotes matrix trace, represent the expectation and covariance with respect to x.

[00138] Lemma 2: Canonical form of a multivariate Gaussian distribution

[00139] The canonical representation of a multivariate Gaussian (Normal) distribution of can be derived as denotes matrix determinant

[00140] Lemma 3: Product of two multivariate Gaussian distributions

[00141] If are two independent random vectors with Gaussian distributions, i.e., then the product of these two distributions can be obtained as follows using the canonical form presented in Lemma 2.

[00142] Using the following definitions

[00143] Eq. (30) can be written as

[00144] Comparing Eqs. (29) and (32), it can be observed that the product of the two Gaussian distributions will be a scaled Gaussian distribution with the scale factor of exp and the mean vector x and the covariance matrix Σ as

[00145] Lemma 4: Trace of a quadratic form

[00146] Considering vector and matrix the quadratic form which is a 1 x 1 matrix can be also represented as

[00147] where tr(. ) denotes matrix trace.

[00148] The present disclosure has described one or more preferred embodiments, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly stated, are possible and within the scope of the invention.