Calibration System and Method for Calibrating an Industrial System Model using Simulation Failure

Abstract

A calibration system and method for calibrating a model of dynamics of an industrial system is provided. The calibration method includes simulating the model multiple times with different combinations of parameters within an admissible range of values of the parameters to estimate success or failure of the simulation. Training the calibration system, iteratively, until a termination condition is met for defining a likelihood of failure of the simulation of the model and for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model and their corresponding calibration errors. Further, calibrating, when the termination condition is met, the model with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping.

Description

TECHNICAL FIELD

The present disclosure relates generally to a calibration system and method for calibrating industrial system models, and more particularly to a calibration system and a method for calibrating an industrial system model based on simulation failure.

BACKGROUND

Advances in modeling and computation have resulted in high-fidelity digital twin models capable of simulating the dynamics of a wide range of industrial systems. The digital twin is a virtual representation of a system that is updated from real-time data and uses simulation, machine learning, and reasoning to help in decision-making. The digital twin uses predictive models of dynamical systems designed for simulation, where the predictive models are often formulated as black-box modeling software that prohibits a user from accessing the underlying model structure.

The approach of formulating the predictive models as the black-box model may be used for a variety of reasons, including the emergent complexity that arises in software that is developed and refined over a period of years or decades, limited application programming interfaces to simplify interactions with larger software architectures, or the desire to avoid the wide distribution of trade secrets or proprietary information encoded in the model.

Despite the fact that this model structure is inaccessible to the user, its parameters often need to be estimated to calibrate the model to observed data. Existing parameter estimation algorithms applicable to black-box systems are designed under the assumption that the model can be simulated forward in time over any duration of interest for any set of parameters within an admissible search domain.

However, in practice, this assumption is invalidated by many existing simulation-oriented models due to multi-scale dynamics, significant nonlinearities, and numerically stiff behavior. Even with recent solvers, the forward simulation of such models often takes a significant amount of time due to small step-sizes and event-driven behavior in hybrid continuous/discrete systems. It is also common for the shape of the admissible parameter set that results in valid simulations to be complex and non-intuitive. Thus, parameters selected from user-defined search domains often result in simulation failure. In the worst-case, a candidate parameter set may result in a very slow simulation that culminates in failure, wasting significant compute. Such behavior can present significant barriers when implementing calibration methods on the black-box model.

Therefore, there is a need for an efficient system and/or method that can estimate near-optimal parameters without extensive simulations to avoid the expenditure of significant time and resources without a corresponding increase in simulation performance.

SUMMARY

Accordingly, it is an object of some embodiments to provide a calibration system and a calibration method that learns from simulation failures and incorporate this information to increase the probability of selecting sets of parameters that lead to successful simulations while avoiding the expenditure of significant time and resources.

Some embodiments are based on the realization that high-fidelity digital models (also referred to as “models”), capable of simulating the dynamics of a wide range of industrial systems, often require calibration or the estimation of an optimal set of parameters in some goodness-of-fit sense, to reflect observed behavior of the industrial system. While searching over parameter space (i.e., admissible range of the parameter) for one or more optimal values of the parameter is an inevitable part of the calibration process, models are seldom designed to be valid for arbitrarily large parameter spaces.

Some embodiments are based on the realization that application of existing calibration methods often results in repeated model evaluations over parameters that can cause the simulations to be unreasonably slow or fail altogether. In general, the shape of subregions in the parameter space that could result in simulation failure is unknown.

Such challenges are particularly common in building energy models. For example, an energy model may seek to describe the temporal behavior of occupied buildings with their associated closed-loop space conditioning systems. Time constants for the closed-loop space conditioning systems typically vary over many orders of magnitude, from milliseconds to weeks, and do not exhibit a clear separation of these timescales. These closed-loop space conditioning systems also often demonstrate clear hybrid behavior due to the cycling behavior of equipment, discontinuities in the derivatives of vapor-compression cycles, and discontinuities in the solar heat gains at sunrise and sunset. Nonlinear interactions between parameters can also cause poor simulation behavior, as when a change in the equipment model parameters results in a loss of cooling capacity and a corresponding significant increase in the room temperature for a given set of building model parameters. Thus, admissible parameter search space is very complex and difficult to capture except through a large investment of time by an expert user.

To mitigate the above-mentioned difficulties, some embodiments of present disclosure provides a calibration system that uses a failure-robust Bayesian optimization (FR-BO) algorithm that learns failure regions in the parameter space from online simulations and informs a Bayesian optimization (BO) algorithm to avoid the failure regions while searching for optimal parameters in the admissible parameter search space. This results in acceleration of the optimizer's convergence and prevents wastage of time trying to simulate parameters with high failure probabilities. In particular, the FR-BO uses simulated data to ascertain regions in the parameter space where the system is likely to fail. The regions correspond to combination of values of a plurality of parameters in the admissible range of the parameters (also referred to as “parameter space”). Further, some embodiments of the present disclosure provides an acquisition function based on a constrained BO acquisition functions that allow searching for optimizer candidates that not only fit operational data well but also are unlikely to result in simulation failures. The operational data is indicative of measurements of an operation of the industrial system.

In particular, in some embodiments, the calibration system calibrate a model of dynamics of the industrial system. The model explains changes of a state of the industrial system in response to controlling the industrial system according to a sequence of control inputs to fit the operational data indicative of measurements of the operation of the industrial system including values of the control inputs and corresponding values of the state of the industrial system. The calibration system simulates at least one operation of the industrial system multiple times, each simulation cycle includes execution of the model of the industrial system having different combinations of values of the parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation.

In some embodiments, the calibration system comprises a failure classifier, where the failure classifier is trained to define a likelihood of failure of the operation of the industrial system for the admissible range of values of parameters of the model using training data including the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation. The calibration system further comprises a probabilistic parameter-to-cost regressor, where the parameter-to-cost regressor is trained for mapping between various combinations of different values of the parameters of the model of the industrial system and their corresponding calibration errors.

In some embodiments, the parameter-to-cost regressor is trained iteratively until a termination condition is met. The probabilistic parameter-to-cost regressor is based on a Bayesian optimization technique. To perform an iteration, a calibration acquisition function of a first two order moments of the calibration errors and a combination of parameters having the maximum likelihood of minimizing a calibration error is identified according to the probabilistic parameter-to-cost regressor. The probabilistic parameter-to-cost regressor defines the first two order moments of the calibration errors. In some embodiments, the at least first two order moments may include a mean of the calibration error and a variance of the calibration error (also referred to as “a confidence range”). For example, for a given a combination of different values of the different parameters of the model of thermal dynamics, the probabilistic parameter-to-cost regressor provides not only a calibration error but also a confidence range around the calibration error. In response, the convergence rate of Bayesian optimization increases and the computational burden of the calibration system reduces.

For example, some embodiments aim to select a combination of the different parameters (also referred as ‘data point’) that need to be queried next. In such a case, querying the combination of the different parameters refers to simulation of a state of the industrial system on the model of the industrial system (such as a thermal state in an environment based on a model of thermal dynamics) with the combination of the different parameters and the received values of the control inputs corresponding to the received values of the state of the industrial system. Some embodiments use an acquisition function of the first two order moments of the calibration error to select the combination of the different parameters to query next.

The acquisition function uses the probabilistic mapping provided by the probabilistic parameter-to-cost regressor to select the combination of the different parameters to query next. In some embodiments, the acquisition function is maximized by the calibration system to select a combination of the different parameters having the largest likelihood of being a global minimum at the probabilistic parameter-to-cost regressor. Therefore, the acquisition function is used as a guidance to select the combination of different parameters to query next. To that end, the calibration system selects the combination of different parameters having the largest likelihood of being a global minimum at the probabilistic parameter-to-cost mapping according to the acquisition function of the first two order moments of the calibration errors.

Further, the calibration system determines a likelihood of failure of the operation of the industrial system for an admissible range of values of parameters from the selected combination of the different parameters of the model. To that end, the calibration system uses training data including values of the selected combinations of the parameters labeled with estimated success or failure of the simulation of the model. Based on the result of the likelihood, the probabilistic parameter-to-cost mapping is updated.

In the next iteration, the calibration system selects a new combination of different parameters having the largest likelihood of being a global minimum at the updated probabilistic parameter-to-cost mapping and estimates a likelihood of failure of the operation of the industrial system for the selected values of the new combination of different parameters. Subsequently, the calibration system again updates the updated probabilistic parameter-to-cost mapping based on the estimated calibration error for the new combination of different parameters, using the Bayesian optimization. Likewise, the probabilistic parameter-to-cost mapping is iteratively computed/updated until a termination condition is met. In an embodiment, the termination condition includes a number of iterations defined by a user.

The calibration system monitors if the termination condition is met. When the termination condition is met, the calibration system outputs an optimal combination of different parameters of the model of industrial system having the largest likelihood of being a global minimum at the probabilistic parameter-to-cost mapping.

The optimal combination of different parameters of the model of the industrial system may be determined in an offline manner (i.e., in advance). Additionally or alternatively, in some implementations, the optimal combination of different parameters of the model of the industrial system may be determined online, i.e., in real-time. For instance, the calibration system may be integrated with a Heating, ventilation, and air conditioning (HVAC) system installed to condition the environment, and the calibration system may output the optimal combination of different parameters during the operation of the HVAC system. Further, in some alternate embodiments, the optimal combination of different parameters of the model of thermal dynamics may be determined offline, and the model of thermal dynamics may be updated online. For instance, the optimal combination of different parameters of the model of thermal dynamics is determined offline and submitted to a controller of the HVAC system. The controller may determine control inputs to the actuator of the HVAC system based on the optimal combination of different parameters. However, when the physical structure of the building and/or the actuators of the HVAC system changes, the calibration system may update the model of thermal dynamics online by estimating a new optimal combination of different parameters in real time according to the change in the physical structure of the building and/or the actuators of the HVAC system. Further, the new optimal combination of different parameters may be submitted to the controller of the HVAC system. Consequently, the controller may determine control inputs to the actuator of the HVAC system based on the new optimal combination of different parameters.

In some embodiments, the Bayesian optimization includes the probabilistic parameter-to-cost regressor using a Gaussian process (GP) for providing the probabilistic mapping, and the acquisition function that exploits the probabilistic mapping provided by the probabilistic parameter-to-cost regressor to direct the querying of a consequent combination of the different parameters of the model of the industrial system.

Some embodiments are based on the recognition that inversion and determinant operations typically used in the GP result in cubic complexity with the number of data points, which implies that using the GP in moderate or high-dimensional parameter spaces is not practical, since finding optimal solutions in such spaces typically requires sampling and evaluating the calibration-cost function a large number of times.

Accordingly, one embodiment of the present disclosure provides a calibration system for calibrating a model of dynamics of an industrial system that explains changes of a state of the industrial system in response to controlling the industrial system according to a sequence of control inputs to fit operational data indicative of measurements of an operation of the industrial system including values of the control inputs and corresponding values of the state of the industrial system. The calibration system comprises: at least one processor; and memory having instructions stored thereon that, when executed by the at least one processor, cause the calibration system to: simulate an operation of the industrial system multiple times, each simulation cycle includes execution of the model of the industrial system having different combinations of values of parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation; train a failure classifier defining a likelihood of failure of the operation of the industrial system for the admissible range of values of parameters of the model using training data including the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation; train a parameter-to-cost regressor for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model of the industrial system and their corresponding calibration errors. The parameter-to-cost regressor is trained iteratively until a termination condition is met, and wherein to perform an iteration the at least one processor is configured to: execute an acquisition function of a first two order moments of the calibration errors configured to identify a combination of parameters having the maximum likelihood of minimizing the calibration errors according to the probabilistic parameter-to-cost mapping; update the identified combination of parameters with a failure-robust acquisition function adjusting the identified parameters according to the failure classifier to reduce the likelihood of the failure of the operation of the industrial system with the model having the updated parameters; simulate the operation of the industrial system controlled by the control inputs in the operational data with the model having the updated parameters to produce simulated states of the industrial system; and update the probabilistic parameter-to-cost mapping based on the updated parameters and calibration errors between the simulated states of the industrial system and corresponding states of the industrial system in the operational data; and calibrate, when the termination condition is met, the model of the industrial system with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping according to the acquisition function.

Accordingly, another embodiment of the present disclosure provides a calibration method for calibrating a model of dynamics of an industrial system that explains changes of a state of the industrial system in response to controlling the industrial system according to a sequence of control inputs to fit operational data indicative of measurements of an operation of the industrial system including values of the control inputs and corresponding values of the state of the industrial system. The calibration method comprises simulating an operation of the industrial system multiple times, each simulation cycle includes execution of the model of the industrial system having different combinations of values of parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation; training a failure classifier defining a likelihood of failure of the operation of the industrial system for the admissible range of values of parameters of the model using training data including the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation; training a parameter-to-cost regressor for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model of the industrial system and their corresponding calibration errors. The parameter-to-cost regressor is trained iteratively until a termination condition is met. Further, to perform an iteration the calibration method comprises: executing an acquisition function of a first two order moments of the calibration errors identifying a combination of parameters having the maximum likelihood of minimizing the calibration errors according to the probabilistic parameter-to-cost mapping; updating the identified combination of parameters with a failure-robust acquisition function adjusting the identified parameters according to the failure classifier to reduce the likelihood of the failure of the operation of the industrial system with the model having the updated parameters; simulating the operation of the industrial system controlled by the control inputs in the operational data with the model having the updated parameters to produce simulated states of the industrial system; and updating the probabilistic parameter-to-cost mapping based on the updated parameters and calibration errors between the simulated states of the industrial system and corresponding states of the industrial system in the operational data; and calibrating, when the termination condition is met, the model of the industrial system with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping according to the acquisition function.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A illustrates a working environment of a calibration system for calibrating a model of an industrial system, according to some embodiments of the present disclosure.

FIG. 1B illustrates modules of the calibration system, according to some embodiments of the present disclosure.

FIG. 1C illustrates a graphical representation of mean and confidence range output from a probabilistic parameter-to-cost regressor module of the calibration system, according to some embodiments of the present disclosure.

FIG. 1D illustrates a graphical representation of an acquisition function executed by the calibration system, according to some embodiments of the present disclosure.

FIG. 2 illustrates a block diagram of the calibration system for calibrating the industrial system, some embodiments of the present disclosure.

FIG. 3 illustrates a first table (Table I) including hyperparameters of a Robertson system for execution of a failure-robust Bayesian optimization (FR-BO) method, according to some embodiments of the present disclosure.

FIG. 4A illustrates graphical representation of comparison of true outputs and estimated outputs of the Robertson system, according to some embodiments of the present disclosure.

FIG. 4B illustrates slices of a parameter space with probabilities of failure estimated by a failure classifier module of the calibration system, according to some embodiments of the present disclosure.

FIG. 5A illustrates a scenario including thermal dynamics of a building and a Heating, ventilation, and air conditioning (HVAC) system to be calibrated by the calibration system, according to an example embodiment of the present disclosure.

FIG. 5B illustrates a second table (Table II) including comparison between best estimates of parameters after predetermined number of iterations during calibration and true values of parameters of the building and the HVAC system, according to an example embodiment of the present disclosure.

FIG. 5C illustrates estimated outputs and true outputs associated with the building, according to an example embodiment of the present disclosure.

FIG. 5D illustrates estimated outputs and true outputs associated with the HVAC system, according to an example embodiment of the present disclosure.

FIG. 6 illustrates regret and number of simulation failure of an integrated building system including the HVAC system and the building, according to an example embodiment of the present disclosure.

FIG. 7 illustrates steps of a calibration method for calibrating a model of the industrial system, according to an example embodiment of the present disclosure.

FIG. 8 illustrates steps of a calibration method for refining admissible parameter space of the simulation model based on the failure region, according to an example embodiment.

DETAILED DESCRIPTION

In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure. It will be apparent, however, to one skilled in the art that the present disclosure may be practiced without these specific details. In other instances, apparatuses and methods are shown in block diagram form only in order to avoid obscuring the present disclosure.

As used in this specification and claims, the terms “for example,” “for instance,” and “such as,” and the verbs “comprising,” “having,” “including,” and their other verb forms, when used in conjunction with a listing of one or more components or other items, are each to be construed as open ended, meaning that that the listing is not to be considered as excluding other, additional components or items. The term “based on” means at least partially based on. Further, it is to be understood that the phraseology and terminology employed herein are for the purpose of the description and should not be regarded as limiting. Any heading utilized within this description is for convenience only and has no legal or limiting effect.

Generally, simulation is a very helpful and valuable work tool. It can be used in industrial field allowing a system's behavior to be learnt and tested. The simulation provides a low cost, secure, and fast analysis tool. It also provides benefits, which can be reached with many different system configurations. With advancement in modeling and computation, high-fidelity digital models capable of simulating the dynamics of a wide range of industrial systems are developed. These models often require calibration, or the estimation of an optimal set of parameters. Practical calibration methods are often designed to estimate near-optimal parameters without extensive simulations to avoid the expenditure of significant time and resources without a corresponding increase in simulation performance.

Accordingly, different methods are used for learning parameters to avoid excessive time and resources wasted in simulation of these models. However, learning the parameters from limited data is challenging. For instance, a Bayesian optimization (BO) method is an effective method for learning the parameters based on limited data in a few-shot manner: that is, with markedly fewer evaluations of the cost function (equivalently, model simulations) than population-based methods. Furthermore, Bayesian optimization inherently balances exploration and exploitation and can incorporate non-convex constraints via modified acquisition functions making Bayesian optimization a powerful and easy-to-use learner for model calibration.

Some embodiments of the present disclosure are based on the realization that the model may not be simulated forward in time over any duration of interest for any set of parameters within an admissible search domain. Therefore, the model simulations can fail to complete in the duration of interest.

Further, some embodiments of the present disclosure are based on the realization that simulation data associated with simulation failure of the model can be used to learn optimal values of parameters, in the entire admissible space of the parameters, that have less likelihood of failure. Thus, the model can be trained to select values of a parameter that decreases chances of the simulation failure and consequently decreases the excessive expenditure of resources on the simulation failure.

Accordingly, the present disclosure provides a calibration system based on a failure-robust Bayesian optimization (FR-BO) method. The FR-BO method comprises different modules that use simulated data to ascertain regions in the parameter space where the industrial system is likely to fail. Subsequently, a specific acquisition function is designed based on constrained BO acquisition functions allowing search for optimizer candidates that not only fit the operational data well, but also are unlikely to result in simulation failures. The calibration system with FR-BO is described further with reference to FIG. 1A and FIG. 1B.

FIG. 1A illustrates a working environment 100 of a calibration system 101 for calibrating a model of an industrial system 103, according to some embodiments of the present disclosure. The model of the industrial system 103 mimics dynamics of the industrial system 103. In particular, the model of the industrial system 103 corresponds to a software model (also called as a mathematical model). In the working environment 100 of the calibration system, the calibration system 101 and the industrial system 103 are coupled to each other via a network 105. Components described in the working environment 100 may be further broken down into more than one component and/or combined in any suitable arrangement. Further, one or more components may be rearranged, changed, added, and/or removed.

In some embodiments, the calibration system 100 includes a transceiver 101a configured to exchange information over the network 105. For example, the transceiver 101a receives a model of the industrial system 103 and data indicative of measurements of operations of the industrial system 103, via the network 105. The industrial system 103 may correspond to any system to be controlled such as a heating, ventilating, and air conditioning (HVAC) system.

The industrial system 103 such as the HVAC system includes actuators including an indoor fan, an outdoor fan, an expansion valve actuator, and the like. The actuators may be controlled according to corresponding control inputs, e.g., a speed of the indoor fan, a speed of the outdoor fan, a position of the expansion valve, a speed of a compressor, and the like. Additionally or alternatively, in some implementations, the control inputs may include a value of temperature and/or a value of humidity. In response to controlling the actuators of the HVAC system according to the corresponding control inputs, thermal state in an environment change.

According to another embodiment, the control inputs are determined based on different parameters of the model of the industrial system 103. The model of the industrial system 103 is a digital twin model of dynamics of the operation of the industrial system 103, such that the model of the industrial system 103 explains the change of the state of the industrial system 103. Thus, the control inputs are variable, and the values of the parameters are fixed. For example, for a digital twin model corresponding to an air-conditioning (AC) system control inputs may comprise temperature required in a room (e.g., 23 degrees), speed of the fan, and the likes which are variable. However, a parameter may comprise a volume of the room to be air conditioned, where the volume of the room is fixed.

For example, in the HVAC system, the different parameters of the model of thermal dynamics define a physical structure of one or combination of the building, the actuators of the HVAC system, and an arrangement of the HVAC system to condition the environment. For instance, the different parameters of the model of thermal dynamics include parameters of a building, using the HVAC system, such as a thickness of a floor of the building, an infrared emissivity of a roof of the building, a solar emissivity of the roof of the building, an airflow infiltration rate, interior room air heat transfer coefficient (HTC), exterior air HTC, and the like. Additionally, the different parameters of the model of thermal dynamics may include HVAC parameters such as an outdoor HEX heat transfer coefficient (HTC) adjustment factor, an indoor HEX HTC adjustment factor, an indoor HEX Lewis number, an outdoor HEX vapor HTC, an indoor HEX vapor HTC, an outdoor HEX liquid HTC, an indoor HEX liquid, an outdoor HEX 2-phase HTC, an indoor HEX 2-phase HTC, and the like.

In some embodiments, the network 105 may be a wireless communication network, such as cellular, Wi-Fi, internet, local area networks, or the like. In some alternative embodiments, the network 105 may be a wired communication network. The calibration system 101 may execute the model of the industrial system 103 multiple times to obtain simulation data. The calibration system 101 uses the FR-BO method to learn, from simulation data, values of parameters for different combinations that increases a likelihood of failure of simulation of the model. In some embodiments, the calibration system 101 may be pre-trained on the simulation data of the model to determine optimal combination of values of the different parameters that can lead to successful simulation.

Further, the transceiver 101a transmits the optimal combination of the different parameters to the industrial system 103, via the network 105. Additionally, in some embodiments, the industrial system 103 may include a controller that determines control inputs for actuators of the industrial system 103, based on the optimal combination of the different parameters received from the calibration system 101. Further, states of the actuators of the industrial system 103 may be controlled based on the control inputs in order to obtain output or a specific state of the industrial system 103 desired by a user.

Further, a mathematical formulation used for modelling the industrial system 103 and challenges faced during determination of optimal values of parameters for combinations of different parameters of the model are elaborated as follows.

Assume that a general model of a dynamical system (also referred to as “the industrial system 103”) is denoted by equation (1):

y _0:T =M _T(θ)   (1)

where y_0:T represents output of a model M _T(θ) of the dynamical system, the constant parameters of the model are described by θ∈Θ⊂^nθ. Further, assume that the admissible set of parameters Θ is known. For instance, Θ could denote a set of upper and lower bounds on parameters obtained from archived data or domain knowledge. The model M_T(θ) is a like a black-box, where a user may not be able to tune parameters of the model M_T(θ). Therefore, range of Θ of the parameters may be purely a guess. Consequently, the range of Θ is not tight around the true parameter set.

Further, the output vector) y_0:T∈^ny×T comprises all measured outputs from the dynamical system obtained over a time period [0, T]. The model M_T(θ) is simulated forward with a fixed (and admissible) set of parameters θ that yields a vector of outputs y_0:T:=[y₀ y₁ . . . y_t . . . y_T ], with each output measurements y_t∈^ny.

In an example embodiment, a model of a state-space representation of a nonlinear dynamical system with state x, initial condition χ₀, control input u, and parameters θ₁ and θ₂ can be modeled as:

{dot over (x)}=ƒ_ODE(χ, u, θ₁),{dot over (y)}=h_ODE(χ, u, θ₂)

Given a set of parameters θ:={θ₁}∪{θ₂}, the nonlinear dynamical system of ordinary differential equations is numerically integrated (i.e., simulated) forward for t∈[0, T], and subsequently the sequence of outputs y_0:T is obtained. Thus, the system dynamics of the dynamical systems (linear and/or nonlinear) can be represented by a parameter-to-output map as in equation (1).

In some scenarios, a class of dynamical systems is modeled by using differential algebraic equations:

0=ƒ_DAE({dot over (x)},x,u,θ₁), {dot over (y)}=h_DAE(x,u,θ₂)

can also be modeled by equation (1).

It is an objective of some embodiments to estimate parameters θ* such that the modeling error (ε) represented mathematically as in equation (2) below can be minimized:

εy*_0:T−M_T(θ*)   (2)

where y*_0:T denotes the measured outputs collected from a real system, and M_T(θ*) denotes the estimated outputs from the model M_T(θ) using the estimated parameters θ*.

The parameters that minimize the modeling error is estimated by optimizing a calibration cost function Ky*_0:T, M_T(θ)) to yield the optimal parameters

$\begin{array}{cc}{\theta }^{*}=\underset{\theta \in \Theta }{\mathrm{argmin}}J\left(y_{0:T}^{*},M_T\left(\theta \right)\right)& \left(3\right)\end{array}$

To that end, some embodiments use BO method which is effective at finding global optima of functions of which gradients are not available and are expensive to evaluate, as is the case in black-box model calibration.

Some embodiments are based on the realization that conventional methods that are used for solving equation (3) assume that M_T(θ) exists for every θ∈Θ, which implies that the model M_T(θ) can be simulated from the time-span of interest [0, T] for any parameter in the admissible parameter space Θ. However, this is not always the case and simulations of the model M_T(θ) can fail to complete in a timespan of interest.

Simulation failure includes at least one case of: (i) a numerical integration schemes terminates prematurely due to parameter-dependent stiffness in the underlying dynamics; (ii) the underlying dynamics do not have a solution due to parameters not adhering to basic validity assumptions, for instance, if the underlying system has log(1−θ²) terms and {|θ|>1}⊂Θ; (iii) a subset of the Θ renders the underlying dynamics unstable (e.g., {dot over (x)}=θx² for θ>0) or a controller/estimator designed using an approximation of the model (such as a linearization) makes the closed-loop dynamics unstable; and (iv) the simulation takes exorbitantly long for some parameters so the code is terminated based on heuristics after a prefixed termination time.

In some cases, models tend to exhibit simulation failures in a failure region Θ_F⊂Θ and the failure does not always occur instantaneously. For instance, in the case (iv), the failure is flagged after a designated termination time that could be large. Consequently, data-driven algorithms that have been designed agnostic to simulation failure could potentially continue to compute optimizer candidates that reside in the failure region Θ_F. In such cases, the algorithms may be deteriorated in performance and lead to large amounts of computational resources and CPU time being wasted.

To overcome the above-mentioned problems, the present disclosure provides the calibration system 101 that is based on a data-driven parameter estimation framework. The calibration system 101 learns from simulation failures and incorporates information from the simulation failure(s) to increase the probability of selecting sets of parameters that lead to successful simulations, thereby optimizing equation (3) without wasting computational resources and CPU time. The data-driven framework is based on the critical assumption that estimated parameters (θ*) do not belong to the failure region (Θ_F) (i.e., θ*∉Θ_F).

The calibration system 101 uses the failure-robust Bayesian optimization approach (implemented using the FR-BO module 107), where initially a probabilistic classifier is designed that estimates the failure region Θ_F from simulation data obtained by sampling within Θ. Further, an entropy-based active learning method is used to reduce sample complexity for evaluation of different functions (such as ƒ_ODE(x, u, θ₁), h_ODE(x, u, θ₁), or the like) to estimate the failure region Θ_F. After estimating {circumflex over (Θ)}_F and the failure region Θ_F, the probabilistic classifier provides probabilities of simulation failure over the entire parameter space of interest Θ. These probabilities are embedded into a Bayesian optimization framework using a constrained acquisition function, ensuring that optimizer candidates are biased to reside outside the estimated failure region {circumflex over (Θ)}_F. Thus, the probabilistic classifier is trained to suggest optimizer candidates that lie within the set difference Θ\ΘF with high probability.

FIG. 1B illustrates modules of the calibration system 101, according to some embodiments of the present disclosure. The calibration system 101 is used for calibrating the model of dynamics of the industrial system 103. The model is calibrated in such a way that the model can explain changes of a state of the industrial system 103 in response to controlling the industrial system 103 according to a sequence of control inputs to fit operational data. The operation data is indicative of measurements of an operation of the industrial system 103 including values of the control inputs and corresponding values of the state of the industrial system 103.

As shown in the FIG. 1B, the calibration system 101 comprises a simulation model (M_T(θ)) 115 which is used to simulate an operation of the industrial system 103 multiple times. In some cases, the simulation model 115 may be provided to the calibration system 101 by a user. Further, the simulation model 115 is simulated multiple times for different combinations of parameters, where values of parameters may be different during each simulation cycle. The simulation model obtains values θ_j of different combination of parameters for each simulation cycle from a database 113, where the database 113 is configured to store operational data D_j comprising different values of each parameter in their admissible range. Each simulation cycle includes execution of the simulation model 115 with different combinations of values θ_j of the parameters selected within the admissible range of values of the parameters. In an example embodiment, the values θ_j for simulating the simulation model 115 may be provided by the user.

Results of the simulations comprises at least one of successful simulation of the simulation model 115 for the selected values of the combination of the parameters or failure of the simulation of the simulation model 115 for the selected values of the combination of the parameters. Further, the results of the simulations of the simulation model 115 are provided to a simulation failure label module 111 of the calibration system 101. The simulation failure label module 111 is configured to label () different combinations of parameters used in multiple simulations of the simulation model 115 based on success or failure of the simulation. For example, the simulation failure label module 111 may be configured to label a first set of combinations of parameters with corresponding values as successful simulation when the simulation is successful and a second set of combinations of parameters with corresponding values as unsuccessful simulation when the simulation is unsuccessful.

Further, based on the result of the calibration, a calibration error is calculated. The calibration system 101 is configured to minimize calibration error to obtain values for optimal combination of parameters. To that end, the calibration system 101 comprises a calibration cost function module 119 configured to calculate the calibration error (J_j) (also referred to as “cost function”) between the simulated states of the model (M_T(θ_j)) of the industrial system 103 and corresponding states of the industrial system 103 in the operational data i.e., true system output data (y*_0:T) 117.

The calculated cost function (J_j) and the labels () generated by the calibration cost function module 119 and the simulation failure label module 111, respectively, are provided to the FR-BO module 107. The FR-BO module 107 is trained iteratively based on the cost function (J_j) and the labels () to obtain a probabilistic mapping between various combinations of different values of the parameters of the model and their corresponding cost functions. The FR-BO module 107 is executed iteratively until a termination condition is met, where the termination condition may be pre-defined by the user. Further, during each iteration, the calibration system 101 determines optimal values of combinations of different parameters that lead to successful simulations, where the optimal values of combinations of different parameters have the largest likelihood of minimizing the calibration error at the probabilistic mapping. Thus, FR-BO module 107 avoids the simulation failure of the model (M_T(θ)) by determining the optimal values of combinations of different parameters.

To obtain the probabilistic parameter-to-cost mapping, the FR-BO module 107 uses a parameter-to-cost regressor module 107c, where the parameter-to-cost regressor module 107c is trained iteratively, until the termination condition is met, to map between various combinations of different values of the parameters of the model (M_T(θ_j)) and corresponding calibration errors. The parameter-to-cost regressor module 107c is further configured to determine a mean μ (a first order moment) and a standard deviation σ (a second order moment) of the calibration errors, which is described further with reference to FIG. 1C.

Further, based on the labels , the calibration system 101 is configured to estimate likelihoods of simulation success and failure on Θ of different values of the parameters by learning the set Θ_F. To that end, the calibration system 101 comprises a failure classifier module 107b, where the failure classifier module 107b is trained, in each iteration “j” of the FR-BO, to define a likelihood of failure (also referred to as “probability of failure PF”) of the operation of the industrial system 103 for the admissible range Θ of values (θ_j) of parameters of the simulation model 111 using training data. The training data includes the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation. Thus, based on the different combination of the parameters, the failure classifier module 107b estimates failure region Θ_F. The failure region Θ_F is a subregion of the admissible parameter space Θ where the simulation model 111 will fail to simulate for one or more combinations of parameters.

Based on the estimation of the failure region Θ_F, the simulation model 111 is refined by refining the admissible range (also referred to as “admissible parameter space”) Θ of the parameters.

In some embodiments, the calibration system 101 obtains predetermined conditions for the failure region Θ_F corresponding to the simulation model 111. The calibration system 101 is further configured to determine whether the failure region Θ_F of the refined simulation model 111 satisfies the predetermined conditions for the failure region. The predetermined conditions for the failure region Θ_F may define size of the failure region Θ_F within the admissible parameter space Θ corresponding to the different combinations of the parameters. When it is determined that the failure region of the refined model does not satisfy the predetermined condition, for example, the failure region is not equal to or less than the predefined size of the failure region Θ_F, the calibration system 101 iteratively refines the model 111 and re-estimate the failure region until the estimated failure region Θ_F satisfies the predetermined condition for the failure region Θ_F.

The calibration system 101 further comprises an acquisition function module 107d that uses an acquisition function to choose parameters such that the failure region Θ_F is avoided. To that end, the acquisition function module 107d obtains the first two order moments of the calibration errors from the parameter-to-cost regressor module 107c and executes the acquisition function for the first two order moments. The acquisition function is executed to identify a combination of parameters having the maximum likelihood of minimizing a calibration error according to the probabilistic parameter-to-cost mapping. In some embodiments, the acquisition function of the first two order moments of the calibration error is used to select the combination of the different parameters to query the next combination of the different parameters.

In particular, the acquisition function module 107d uses the probabilistic mapping provided by the parameter-to-cost regressor module 107c to select the combination of the different parameters to query next, which is described further with reference to FIG. 1D.

In an embodiment, the acquisition function is maximized, in each iteration “j” of the FR-BO, by the calibration system 101 to select a combination of the different parameters having the largest likelihood of being a global minimum at the parameter-to-cost mapping provided by the parameter-to-cost regressor module 107c. Therefore, the acquisition function is used as a guidance to select the combination of different parameters to query next. To that end, the calibration system 101 selects the combination of different parameters having the largest likelihood of being a global minimum at the probabilistic parameter-to-cost mapping according to the acquisition function of the first two order moments of the calibration errors. The acquisition function module 107d uses at least one of the acquisition functions: expected improvement (EI), the lower confidence bound, or entropy search. In a preferred embodiment, the acquisition function module 107d uses the EI acquisition function.

Further, the calibration system 101 comprises a failure-robust acquisition function module 107a that is configured to adjust values of the parameters in the identified combination of parameters based on the output of the acquisition function module 107d and output of the failure classifier module 107b. Thus, the failure-robust acquisition function module 107a incorporates probabilistic output into the acquisition function to yield failure-robust EI (FREI) acquisition function. Further, the updated values of the combination of parameters are used for training the FR-BO module 107 in next iteration. The updated values of the combination of the parameters (θ_j) are stored in the databases 113 and are used for simulating the simulation model 115.

In this way the FR-BO module 107 is trained iteratively on the simulation data until the termination condition is met.

Failure Robust Bayesian Optimization (FR-BO) Algorithm

The FR-BO module 107 implements the FR-BO algorithm using modules 107a -107d (as shown in FIG. 1B). The FR-BO module 107 learns the failure region from data obtained during simulations via a probabilistic classifier such as a scalable variational Gaussian process classifier (VGPC). The probabilities of simulation failure or simulation success are used via active learning method to accelerate the Bayesian optimization step, and further, to guide where best to simulate the dynamical model M_T(θ) to get better estimates of the failure region boundaries. Further, information from the VGPC is incorporated via a constraint weighted acquisition function to ensure that parameters are chosen avoiding failure regions.

To that end, the FR-BO module 107 is trained to learn the failure regions. Since the admissible parameter search domain Θ is known, by sampling on admissible parameter space (or Θ ), the training data is obtained for learning the failure region subset Θ_F. Further, most informative (in an information-theoretic sense) parameters θ are selected from the admissible parameter search domain Θ (i.e., θ∈Θ). For instance, the training data is initially available obtained by random sampling on Θ. Based on the training data, the FR-BO module 107 is trained iteratively until the termination condition is met i.e., the failure classifier 107b and the parameter-to-cost regressor 107c are trained iteratively.

At the j-th iteration of training the failure classifier module 107b, the training data D_j comprises a sequence of parameters θ_[0:j], a sequence of corresponding simulation failure labels _[0:j], and cost function values J_[0:j].

D _j=θ_[0:j]×_{[0: j]}×J_[0:j]

Each failure label is denoted as “+1” if there is simulation failure and donated as “−1” if no simulation failure. For failed simulations, the corresponding cost function value is set to an undefined value (or null value), e.g., NaN. If the simulation is successful, the cost function yields a real-valued scalar.

Further, at the j-th iteration of learning the failure region, θ and labels of the dataset D_j are utilized to set a Gaussian process prior (i.e., prior probability distribution function) at the observed parameter sets. This can be written as ϕ˜N(0, K_j) , where ϕ_j is the prior function using D_j and with a user-specified kernel function K(⋅,) such as a squared exponential kernel or a Matern Kernel.

$\begin{array}{cc}{K}_j=\left[\begin{array}{ccc}\left({\theta }_0,{\theta }_0\right)& \dots & \left({\theta }_0,{\theta }_j\right)\\ ⋮& \ddots & ⋮\\ \left({\theta }_j,{\theta }_0\right)& \dots & \left({\theta }_j,{\theta }_j\right)\end{array}\right]& \left(4\right)\end{array}$

To perform classification of the combination of parameters into simulation success and simulation failure with this prior, the function ϕ is transformed through a squashing function such as the cumulative density function of a zero-mean unit-variance normal distribution γ(⋅):=N(⋅|0, 1), given by Γ(z)=∫_−∞^zγ(α)dα. Consequently, a Bernoulli distribution can be used to represent a likelihood function conditioned on the transformed data as follows:

(_j|Γ(ϕ_j))=Γ(ϕ_j)·(1−Γ(ϕ_j)

The joint distribution of and ϕ thus becomes

p(,ϕ)=Π_r=1^j(_j|Γ(ϕ_j)) N(0, K_j)   (5)

To optimize hyperparameters and perform inference two more distributions are required: the marginal likelihood PF(, j) and the posterior p(ϕ|, j). Both of these distributions require the inversion of the j×j kernel matrix (represented in equation (4)). However, the inverted kernel matrix incurs cubic complexity and does not scale well to the large values of j that may be required for FR-BO to compute good solutions.

To that end, approximation methods that leverage pseudo-inputs are used. The approximation methods are also known as inducing points. Inducing points {tilde over (θ)}∈Θ are design variables that are augmented with the latent variables ϕ_j that also respect the Gaussian prior and yields a joint distribution

$\left(\varphi ,\tilde{\varphi }\right)\sim N\left(0,\left[\begin{array}{cc}{K}_j& {\tilde{K}}_{jm}\\ K_{jm}^T& {\tilde{K}}_m\end{array}\right]\sim \right)$

{tilde over (K)}_jm denotes the covariance matrix computed by evaluating the kernel across j data points and m inducing inputs, while {tilde over (K)}_m denotes the covariance matrix computed by evaluating the kernel on all pairs of the inducing inputs. By exploiting the properties of the Gaussian distribution, the joint distribution of the latent variables and the inducing variables are rewritten as:

p(,ϕ,{tilde over (ϕ)})=p(|ϕ)p(ϕ|{tilde over (ϕ)})p({tilde over (ϕ)})

Further, to get a variational approximation of the likelihood, the following inequality is used:

log p(|{tilde over (ϕ)})≥E_{p( ϕ|{tilde over (ϕ)})}[log p(|ϕ)]

Further, to get the variational lower bound, a variational distribution q is defined

log p()≥E_q(ϕ)[log p(|ϕ)]−KL[q({tilde over (ϕ)})∥p({tilde over (ϕ)})]

The optimal hyperparameters for the VGPC can be obtained by minimizing the loss function formed by the negative of the right hand side of inequality in equation (6) by using quadrature methods. If it is assumed that q˜N({tilde over (ϕ)}|{tilde over (μ)}, {tilde over (Σ)}), then

q(ϕ)=N(L{tilde over (μ)},K_j+L({tilde over (Σ)}+{tilde over (K)}_m)),   (7)

with L={tilde over (K)}_jm{tilde over (K)}_m⁻¹, which is an m×m matrix, and eventually m<<j, so this matrix is cheaper to invert, which makes this method scalable.

In some embodiments, for inference at a set of test points {θ_*}, the set of test points {θ_*} is transformed into {Θ_*}, and the approximate posterior is then given by p(Θ_*|)=∫p({tilde over (Θ)}_*|{tilde over (Θ)})q({tilde over (Θ)})d{tilde over (Θ)}, which can be computed in a manner similar to equation (7). Initially, when the failure region is not correctly estimated, informative elements of {θ*} are selected that can yield good estimates of Θ_F without exhaustive sampling. To this end, the calibration system 101 uses an active learning strategy, where the most informative θ*_j∈θ_* is selected based on the maximum entropy of the posterior distribution. Since the mean and variance of the posterior p(Θ_*|) is computed for each parameter in θ_*, the entropy is computed (assuming q is Gaussian) and fix

${\theta }_j^{*}=\underset{{\theta }^{\prime }\in {\theta }_{*}}{\mathrm{argmax}}\frac{1}{2}\log \left(2\pi \mathrm{Var}\left({\theta }^{\prime }\right)\right)$

θ*_j is then evaluated by simulating M_T(θ*_j) to ascertain _j+1 and J_j+1, which yields the updated dataset D_j+1 and the process iterates till a termination criterion such as a maximum number of iterations is achieved.

Thus, the active learning method accelerates the Bayesian optimization step, and guides the calibration system 101 where best to simulate the dynamical model M_T(θ) to get better estimates of the failure region boundaries.

Classical BO methods assume the presence of one global optimum, and smoothness of the θ to J map. Since J is typically assumed to be continuous, the data at the j-th iteration is leveraged to construct a surrogate GP model of the reward, given by

Ĵ _j :=GP(μ(θ;D_j),σ(θ,θ′;D_j))   (8)

where μ(⋅) is the predictive mean function, and σ(⋅,) is the predictive variance function, and D_j containing {θ_[0:j], J_[0:j]} is the dataset collected thus far. Typically, the variance is expressed through the use of kernel functions.

At the j-th learning iteration, for a new query sample θ∈Θ, the GP model predicts the mean and variance of the reward to be μ(θ)=k_j(θK_j⁻¹J_0:j and σ(θ)=(θ, θ)−k_j(θ) K_j⁻¹k_j(θJ_0:j, where k_j(θ)=[(θ₀, θ)(θ₁, θ) . . . (θ_j, θ)], and K_j is defined in equation (4)

The accuracy of predicted mean and variance is strongly linked to the selection of the kernel and the best set of hyperparameters such as length scales and variance parameters of the kernels and estimated noise. These hyperparameters are obtained by maximizing the log marginal likelihood function (MLL)

$-\frac{1}{2}\log ❘K_j❘-\frac{1}{2}{J\left(\theta \right)}^T{K}_j^{-1}J\left(\theta \right)-\frac{n_{\theta }}{2}\log 2\pi$

In some embodiments, inducing variables as in the VGPC can be used to significantly improve scalability of the GP.

In Bayesian optimization, the mean and variance of a surrogate model Ĵ_j (also referred to as “surrogate GP model”) in equation (8) to construct an acquisition function to inform the selection of a θ_j that increases the likelihood of minimizing the current best cost. To this end, the incumbent

${\hat{J}}_j^{*}:=\underset{\theta \in \Theta }{\min }\mu \left(\theta ;D_j\right)$

is computed and used to define an expected improvement (EI) acquisition function that has the form

$\mathrm{EI}\left(\theta ,j\right)=\{\begin{array}{c}\sigma \left(\theta \right)\gamma \left(z\right)+\left({\hat{J}}_j^{*}-\mu \left(\theta \right)\right)\Gamma \left(z\right)\mathrm{if}\sigma \left(\theta \right)>0,\\ 0\end{array}$

where

$z=\frac{{\hat{J}}_j^{*}-\mu \left(\theta \right)}{\sigma \left(\theta \right)},$

and γ(.), Γ(.) are the probability distribution function (PDF) and the cumulative distribution function (CDF) of the zero-mean unit-variance normal distribution, respectively. Thus, the Bayesian optimization includes the acquisition function that exploits the probabilistic mapping provided by the probabilistic GP surrogate model to direct querying of consequent combination of the different parameters.

In the j-th iteration of learning, the data D_j is used to construct the EI acquisition function using the surrogate Ĵ_j. Subsequently, the optimizer candidate is computed

$\begin{array}{cc}{\theta }_{j+1}=\underset{\theta \in \Theta }{\mathrm{argmax}}\mathrm{EI}\left(\theta ,j\right)& \left(9\right)\end{array}$

which serves as the parameter estimate θ in equation (1) in the next BO iteration. In some cases, other acquisition functions such as the lower confidence bound or entropy search can be used instead of EI.

Further, since the VGPC generates a probabilistic output associated with failure and success of simulation, the probabilistic output can be directly incorporated into the acquisition function. This yields the failure-robust EI (FREI) acquisition function:

FREI(θ,j)·(1−PF(θ,j))   (10)

where PF(θ,j) is the likelihood of failure calculated by training the VGPC using data up to the j-th iteration, and then evaluating the likelihood of the VGPC at θ. If the VGPC algorithm does not find any θ such that P(θ)>0, then the FREI acquisition function (10) is zero for every θ∈Θ and future candidates are selected randomly until at least one 0 is found which allows for a successful simulation.

To that end, both successful simulations and failure simulations are observed, and the VGPC is trained on a non-trivial classification problem. In such a case, the higher the value of P(θ, j), the higher the probability that a particular candidate is selected, so long as its expected improvement is high as well. The multiplicative nature of the components in (10) seeks to ensure that neither one component can outweigh the other, and candidates will be selected only if they are both a candidate for optimization and feasibility (i.e., is expected to yield a successful simulation). Along the same lines as (9), the FR-BO selects the next optimizer candidate as follows:

$\begin{array}{cc}{\theta }_{j+1}=\underset{\theta \in \Theta }{\mathrm{argmax}}\mathrm{FREI}\left(\theta ,j\right)& \left(11\right)\end{array}$

Thus, in the FREI, likelihood is estimated off the binary labels , in contrast with constrained BO, where probabilities are estimated based on continuous slack variables obtained from the constraints. FREI uses Gaussian process (GP) proxies for constraint modeling for a classifier that classifies combinations of parameters into at least one of: parameters causing simulation failure or parameters causing simulation success.

Some embodiments are based on the realization that retraining the VGPC is expensive after the collection of each new data sample. However, as long as the VGPC has been trained initially with some data, VGPC can be retrained infrequently. In such a case, frequency of retraining the VGPC is problem dependent. In some embodiments, the VGPC is retrained when the FR-BO operation is stuck at a local optimum for a pre-decided number of iterations.

FIG. 1C illustrates a graphical representation of mean and confidence range output from a probabilistic parameter-to-cost regressor module 107c of the calibration system, according to some embodiments of the present disclosure. In particular, FIG. 1C illustrates a graph of a parameter (x) on an X-axis and a dependent function (f(x)) on a Y-axis. The graph indicates changes in the dependent function f(x) for different values the parameter X within the admissible range of the parameter X (in this case [−1 to 3]). For the ease of explanation, the graph in FIG. 1C includes only one parameter and only one dependent function. However, in real-life models of the industrial system 103, the graph will comprise a plurality of parameters and a plurality of dependent function.

Referring to FIG. 1C, dots (such as dots 21) within the graph indicates samples/observations, a curve 123 indicates the mean, and shaded area 125 indicates the confidence range. Thus, the graph provides information about estimated value of the dependent function f(x) for an estimated value of the parameter X. The calibration system 101 uses the graph to determine one or more failure regions within the admissible range of the parameter X that may result into failure of the simulation. From the graph of FIG. 1C, it can be determined that there are high chances of the failure of the simulation by selecting values of X in the range 1.25 onwards. This region within the admissible space (or range) of the X is referred to as the failure region because it is high likely that by selecting values of X within the failure region (i.e., from 1.25 to 2) the simulation of the model may fail.

In some embodiments, based on a graph (like the graph in FIG. 1C) corresponding to a model with a plurality of parameters, the calibration system 101 may determine the failure regions as a combination of values of the plurality of parameters that increases the likelihood of failure of the operation of the industrial system 103. The calibration system 101 uses the FR-BO module 107 to learn failure regions of different parameters used in the model of the industrial system 103 and thus, accurately estimates values of the different parameters that lead to successful simulation. In this way, the calibration system 101 saves excessive time and other resources on the simulation.

FIG. 1D illustrates a graphical representation of an acquisition function executed by the calibration system 101, according to some embodiments of the present disclosure. In particular, in FIG. 1D, there is shown a graph including a curve 127 of the acquisition function for selecting the combination of the different parameters which are to be queried next. The acquisition function uses the probabilistic mapping provided by the probabilistic parameter-to-cost module 107c to select the combination of the different parameters to be queried next. In some embodiments, the acquisition function is maximized by the calibration system 101 to select a combination of the different parameters having the largest likelihood of being a global minimum 129 at the probabilistic mapping, for querying. Therefore, the acquisition function is used as a guidance to select the combination of different parameters to query next.

FIG. 2 illustrates a block diagram of the calibration system 200 for calibrating the industrial system 103, according some embodiments of the present disclosure. The calibration system 200 can have a number of interfaces connecting the calibration system 200 with other systems and devices. For example, a network interface controller (NIC) 201 is adapted to connect the calibration system 200, through a bus 203, to a network 205. Through the network 205, either wirelessly or through wires, the calibration system 200 may receive the data 207 indicative of measurements of the operation of the industrial system including values of the control inputs to the actuators of the industrial system and values of a state of the industrial system caused by the operation of the industrial system according to the values of the control inputs.

For example, the industrial system 103 corresponds to an HVAC system. The calibration system 200 may wirelessly receive, via the network 205, values of control inputs to actuators of the HVAC system and values of thermal state at locations of an environment caused by the operation of the HVAC system according to the values of the control inputs. In some cases, the data 207 indicative of the measurements of the operation of the HVAC system and the values of the thermal state at locations of the environment may be received via an input interface 209.

The calibration system 200 includes a processor 211 configured to execute stored instructions, as well as a memory 213 that stores instructions that are executable by the processor 211. The processor 211 can be a single core processor, a multi-core processor, a computing cluster, or any number of other configurations. The memory 213 can include random access memory (RAM), read only memory (ROM), flash memory, or any other suitable memory systems. The processor 211 is connected through the bus 203 to one or more input and output devices. Further, the calibration system 200 includes a storage device 215 adapted to store different modules storing executable instructions for the processor 211. The storage device 215 can be implemented using a hard drive, an optical drive, a thumb drive, an array of drives, or any combinations thereof.

The storage device 215 is configured to store the FR-BO module 107, where the FR-BO module 107 is configured to determine optimal values of combinations of different parameters to avoid the simulation failure of simulation model (not shown in FIG. 2) corresponding to the industrial system 103. To that end, the FR-BO module 107 is trained based on the training data including various combinations of different parameters resulting in simulation failure or simulation success.

The FR-BO module 107 comprises the parameter-to-cost regressor module 107c to obtain the probabilistic parameter-to-cost mapping. The parameter-to-cost regressor module 107c is trained iteratively, until the termination condition is met, to map between various combinations of different values of the parameters of the model and corresponding calibration errors. The parameter-to-cost regressor module 107c is further configured to determine the first two order moments (i.e., the mean μ and the standard deviation σ) of the calibration errors.

The FR-BO module 107 is further configured to estimate likelihoods of simulation success and failure on admissible range (or space) of different values of the parameters by learning the set failure regions. To that end, the FR-BO module 107 uses the failure classifier module 107b, where the failure classifier module 107b is trained to define a likelihood of failure of the operation of the industrial system 103 for the admissible range of values of parameters of the model using the training data. The training data includes the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation.

The FR-BO module 107 further comprises the acquisition function module 107d that uses an acquisition function to select parameters such that the failure region corresponding to the selected parameters is avoided. To that end, the acquisition function module 107d obtains the first two order moments of the calibration errors from the parameter-to-cost regressor module 107c and executes the acquisition function for the first two order moments to identify a combination of parameters having the maximum likelihood of minimizing the calibration error according to the probabilistic parameter-to-cost mapping. In some embodiments, the acquisition function of the first two order moments of the calibration error is used to select the combination of the different parameters to query the next combination of the different parameters.

Further, the FR-BO module 107 comprises the failure-robust acquisition function module 107a that is configured to adjust values of the parameters in the identified combination of parameters based on the output of the acquisition function module 107d and output of the failure classifier module 107b. Thus, the failure-robust acquisition function module 107a incorporates probabilistic output (obtained from the failure classifier module 107b) into the acquisition function to yield failure-robust EI (FREI) acquisition function. Further, the updated values of the combination of parameters are used for training the FR-BO module 107 in next iteration. The updated values of the combination of the parameters are used for simulating the simulation model in the next iteration.

Additionally, the calibration system 200 may include an output interface 217. In some embodiments, the calibration system 200 is further configured to submit, via the output interface 217, the optimal combination of the different parameters of model of the industrial system 103 to a controller 219 of the industrial system 103. The controller 219 is configured to generate control inputs to the actuators of the industrial system 103 based on the optimal combination of the different parameters of the model.

In some embodiments, the calibration system 101 is used to estimate optimal values for different parameters of the Robertson system, where the Robertson system is a stiff nonlinear system of chemical kinetics. The Robertson system is modeled as follows:

{dot over (x)} ₁=−θ₁x₁+θ₂x₂x₃ y=x

{dot over (x)} ₂=θ₁x₂−θ₂x₂x₃−θ₂x₂ ² 0=x₁+x₃+x₃−1

{dot over (x)} ₃=−θ₃x₂ ²,

where x₁, x₂, and x₃ are concentrations of chemical species, and the true parameters of the system are θ₁=4×10⁻², θ₂=10⁻⁴, and θ₃=3×10⁷.

The system of equations is assumed unknown during parameter estimation using FR-BO approach. Ground truth data is collected by simulating the Robertson system forward for T=[0, 20] s, after which T is partitioned into 100 equidistant samples in order to obtain the measured outputs y*_0;T collected from the Robertson system.

Further, the admissible parameter space of parameters associated with the Robertson system is defined in logarithmic space and is given by log₁₀Θ={[−3, −1]×[2, 12]×[5, 15]}z. Thus, θ₁ can lie in the interval [10⁻³, 10⁻¹], θ₂ can lie interval [10⁻², 10¹²], and similarly, θ₃ can lie interval [10⁵, 10¹⁵]. For each θ sampled on Θ, the Robertson system is simulated over T. In an example embodiment, the Robertson system may be simulated using the Radau IIA solver, which is well-suited to the numerical integration of stiff systems.

Further, for each θ, M_T(θ) is computed and a parameter estimation cost function J(y*_0:T, M_T(θ))=log Σ_t=0^T(y*_t−y_tW(y*_t−y_t), where y_t is the output vector obtained from M_T(θ) at the t-th time instant is also estimated. The matrix W is a scaling matrix to ensure that the three output components are of similar magnitudes.

To that end, the scaling matrix W is set to W=diag[1, 5 ×10⁴, 1]. Further, 300 data points are collected using active learning for the VGPC and then run 700 BO iterations, with the VGPC retrained every 100 iterations during BO.

FIG. 3 illustrates a first table (Table I) including hyperparameters of a Robertson system for execution of a failure-robust Bayesian optimization (FR-BO) method, according to some embodiments of the present disclosure. In particular, the first table (Table I) includes hyperparameters used for the VGPC and the GPR learners needed for FR-BO.

After applying FR-BO, the set of parameters that produces the lowest cost is given by {tilde over (θ)}₁=3.717×10⁻², {tilde over (θ)}₂=0.887×10⁴, and {tilde over (θ)}₃=3.03×10⁷, which is satisfactorily close to the true parameters. The effectiveness of the estimated parameters is validated on a time span of [10⁻⁵, 10⁵] to showcase multi-scale dynamics. The similarity of the model outputs and the true outputs is demonstrated in FIG. 4A. Further, the mean squared error between the two trajectories is 9.910×10⁻⁵.

FIG. 4A illustrates graphical representation of comparison of true outputs and estimated outputs of the Robertson system, according to some embodiments of the present disclosure. The comparison is computed in a testing scenario for the time span of [10⁻⁵, 10⁵]. From FIG. 4A, it is clear that the estimated values of the parameters are very close to the real values of the parameters.

Further, two slices of the failure region learned by the VGPC are demonstrated with reference to FIG. 4B.

FIG. 4B illustrates slices of a parameter space with probabilities of failure estimated by a failure classifier module 107b of the calibration system 101, according to some embodiments of the present disclosure. The slice on the left in the FIG. 4B is obtained by fixing θ₂ at its estimated value and allowing θ₁ and θ₃ to vary over the subset of Θ. Similarly, the slice of the right in the FIG. 4B is obtained by fixing θ₁ and allowing the other two parameters θ₂ and θ₃ to vary. Both slices clearly indicate that the VGPC has identified a region where the model simulation is likely to exhibit failure. Further, it is also evident that such a region can have significant curvature. Therefore, even for a low-dimensional dynamical systems, it is advantageous to use advanced nonlinear classifiers like VGPC rather than linear classifiers.

In another example embodiment, the calibration system 101 is used to estimate parameters for dynamical model in a building digital twin. Simulation model of building energy behavior for digital twin application is designed to predict multiscale nonlinear dynamics that result from the temporal interactions between the many different building subsystems.

FIG. 5A illustrates a scenario including thermal dynamics of a building 500 and a Heating, ventilation, and air conditioning (HVAC) 501 system to be calibrated by the calibration system 101, according to an example embodiment of the present disclosure. The HVAC system 501 is installed in the building (also referred to as “building system”) 500 to condition an environment of the building. A model of a vapor compression cycle of the HVAC system 501 comprising a compressor, a condensing heat exchanger, an electronic expansion valve, and an evaporating heat exchanger, is constructed. The vapor compression cycle behavior is dominated by the heat exchangers over time scales of interest. Therefore, the model of a vapor compression cycle uses dynamic models of the heat exchangers and algebraic (static) models of the compressor and the electronic expansion valve. For the ease of describing, a lumped parameter method was used to characterize dynamics of a refrigerant flow in the heat exchangers. In an embodiment, a Helmholtz equation-of-state based model may be used to describe the refrigerant properties, while an ideal gas mixture may be used for the moist air model.

A simple isenthalpic model of the electronic expansion valve and mass flow rate is regularized in a neighborhood of zero flow to prevent a derivative of the mass flow rate from tending toward infinity. A flow coefficient is determined via calibration against experimental data. An operation of the compressor may be described by relating a volumetric efficiency and an isentropic efficiency to a suction pressure, a discharge pressure, and a compressor frequency.

In an example embodiment, building models of the building 500 are constructed from an open-source Modelica Buildings library. A room model is based on a physics-based behavior of fundamental materials and commonly used components, while a zone air model is a mixed air single-node model with one bulk air temperature that interacts with all of the radiative surfaces and thermal loads in the room.

According to an embodiment, the building models may be parameterized by a thickness of a roof 503, a thickness of outside wall 505, and a thickness of floor 507. The thickness of the roof depends on a thickness of concrete layer 509, a thickness of insulation layer 511, and a size of plenum 513. The thickness of floor 507 depends on a thickness of concrete slab 515 and a thickness of the carpet tile 517.

For instance, the building model consists of a one-story residence with nominal 2009 IECC-based construction, where the one-story residence has a floor area, e.g., 112.24 m and is 2.6 m tall, and is oriented along with cardinal directions with a peak occupancy of 3 people. Each outside wall also has a window 519 of size (e.g., 1.52 m ×2.72 m) that admits solar heat gains into spaces of the building 500. A 10 cm thick concrete slab and 2 m of soil below the building 500 is also included to characterize interactions with thermal boundary conditions under the building 500, which was set to a constant 21° C. Additionally, a peaked attic is also included with a maximum height of 1.5 m, so that the building model includes two thermal zones.

Further, the building model is connected to the vapor compression cycle model, and a proportional integral (PI) controller 521 is implemented on a heat pump which used the compressor frequency to regulate a room temperature and the expansion valve position to regulate an evaporator superheat temperature. The PI controller 521 also implemented anti-windup to maintain stability while enforcing minimum and maximum actuator limits. The resulting joint building envelope/HVAC model (i.e., the model of thermal dynamics) is simulated using Atlanta-Hartsfield TMY3 file, and included convective and radiative heat loads of 2 W/m² and a latent load of 0.6 W/m² between hours of 8 AM and 6 PM, with weather-driven disturbances outside of these hours. Further, such a model is exported from Modelica using a Functional Mockup Interface, and resulting functional mockup unit (FMU) is imported into Python using the FMPy package to enable seamless integration of advanced machine learning modules.

Inputs and outputs of the model of thermal dynamics may be chosen to be similar to those which may be observed in a realistic experimental setting. Inputs of the heat pump include a room temperature set-point, an evaporator superheat set-point, and indoor and outdoor fan speeds. Inputs for the building model include convective, radiative, and latent heat loads as well as weather variables provided in TMY3 standard. Such heat loads may be estimated to reasonable accuracy via occupancy detection, load surveys, or other similar methods.

Further, ground-truth data is collected for calibration by simulating the models for a predetermined amount of time. In an example embodiment, the ground truth data is collected by simulating the Modelica model for a period of time (e.g., T=14 days) with 17 parameters of the model set to their true values. The search space of the parameters is ±15% of the true parameters with ±5% random translation of the upper and lower bounds for each parameter. The eight measured output sequences y*_0:Tof the model are collected at 5 minutes intervals, and the FR-BO components have the same hyperparameters as in the Table I (FIG. 3), with two differences: the training iterations are kept at 2000 for the GP, and the kernels are Matern-5/2 for the VGPC.

The GP is initialized by choosing 100 randomly selected parameter samples from within the bounds Θ associated with each parameter. With each of the initial parameters' samples, the model is simulated for the same time interval as the ground truth and obtain the estimated output sequence y*_0:T. Subsequently, the calibration-cost function from equation (3) is evaluated for each of the initial samples with the simulated and measured outputs. The initial collection of parameters and calibration-cost function values is used to construct an initial training set of the GP. In an embodiment, Matern 3/2 kernels with dimension-wise separate length-scales may be used, as the admissible parameter space is not normalized.

Further, the parameters chosen to be representative of the HVAC system 501 and the building 500 are provided in Table II illustrated in FIG. 5B.

FIG. 5B illustrates a second table (Table II) including comparison between best estimates 523 of parameters 525 after predetermined number of iterations during calibration and true values 527 of parameters 525 of the building and the HVAC system, according to an example embodiment of the present disclosure. It can be observed that values of the parameters 525 estimated using FR-BO approach are very close to the true values 527 of the parameters 525.

Further, a comparison of the measured and estimated simulation outputs of the building 500 and HVAC system 501 is shown in FIG. 5C-5D.

FIG. 5C illustrates estimated outputs and true outputs associated with the building 500, according to an example embodiment of the present disclosure. The estimated outputs are provided by the model of the building system 500 calibrated using the FR-BO approach (using the calibration system 101). In the plots of FIG. 5C, the circles are true data points and the continuous lines are the estimated outputs from model simulation. It is clear from the FIG. 5C that the output error for the building system 500 is small, and this is further evident from the second Table II (FIG. 5B), where most parameters have been estimated with high accuracy.

FIG. 5D illustrates estimated outputs and true outputs associated with the HVAC system 501, according to an example embodiment of the present disclosure. The estimated outputs are provided by the model of the HVAC system 501 calibrated using the FR-BO approach (used by the proposed calibration system 101). In the plots of FIG. 5D, the circles are true data points and the continuous lines are the estimated outputs from model simulation. It can be observed from the FIG. 5D that the output error for the HVAC system 501 is small, and this is further evident from the second Table II (FIG. 5B), where most parameters have been estimated with high accuracy.

Further, effectiveness of the proposed calibration system 101 in reducing the amount of time wasted on failed simulations during optimization of an integrated building system (i.e., the building system 500 and the HVAC system 501) is also evaluated.

FIG. 6 illustrates regret and number of simulation failure of an integrated building system including the HVAC system 501 and the building 500, according to an example embodiment of the present disclosure. For the evaluation purpose, the VGPC is allowed to first collect data via active learning for 450 iterations before it begins the optimization process. In the left subplot 601 of FIG. 6, it can be observed that due to the incorporation of the FREI function the number of simulation failures is reduced. This improves the convergence of the Bayesian optimization step, so that the FR-BO converges to its best solution within approximately 700 iterations, which results in a lesser cost.

In the right subplot 603 of FIG. 6, it can be observed that the total time that is used for Bayesian optimization is higher for FR-BO (15 hr vs. 9 hr) and the total time wasted simulating failures is significantly less (15 hr vs. 23 hr), which equates to saving 8 hours of wall time.

FIG. 7 illustrates steps of a calibration method 700 for calibrating a model of the industrial system 103, according to an example embodiment.

At step 701, an operation of the industrial system 103 is simulated multiple times. In each simulation cycle, execution of a model of the industrial system 103 includes different combinations of values of parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation. The industrial system 103 may correspond to any systems that can be modeled such as a building system, an HVAC system, and the like.

At step 703, for defining a likelihood of failure of the operation of the industrial system 103 for the admissible range of values of parameters of the model, the failure classifier 107b is trained. Training data used for training the failure classifier 107b includes the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation.

At step 705, for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model of the industrial system 103 and their corresponding calibration errors, the parameter-to-cost regressor 107c is trained iteratively.

To train the parameter-to-cost regressor 107c iteratively, at step 705a, an acquisition function of a first two order moments of the calibration errors identifying a combination of parameters having the maximum likelihood of minimizing the calibration errors is executed according to the probabilistic parameter-to-cost mapping.

At step 705b, the identified combination of parameters is updated with a failure-robust acquisition function module 107a by adjusting the identified parameters according to the failure classifier module 107b to reduce the likelihood of the failure of the operation of the industrial system 103 with the model having the updated parameters.

At step 705c, the operation of the industrial system103 controlled by the control inputs in the operational data is simulated with the model having the updated parameters to produce simulated states of the industrial system 103. Further, the probabilistic parameter-to-cost mapping is updated based on the updated parameters and the calibration errors between the simulated states of the industrial system 103 and corresponding states of the industrial system 103 in the operational data.

At step 707, the model of the industrial system 103 is calibrated with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping according to the acquisition function.

Some embodiments are based on the realization that the simulation model 111 can be refined based on identification of the failure region. To that end, the calibration system 101 is configured to learn the failure region. The failure region is a subregion of the admissible parameter space where the simulation model 111 will fail to simulate for one or more combinations of parameters. Thus, the failure region is deterministic. To identify the failure regions, the calibration system 101 uses the failure classifier module 107b that is trained to estimate the failure region in the admissible parameter space for different combinations of the parameters. Based on the training, the failure classifier module 107b estimates the failure regions in the admissible parameter space with some likelihood. Thus, the output of the failure classifier module 107b is probabilistic.

Based on the failure regions estimated by the failure classifier module 107b, the calibration system 101 is further configured to refine the admissible parameter space iteratively, until the estimated failure region satisfies predetermined conditions. The predetermined conditions may include size of the failure region. More details regarding the simulation model refinement (or admissible parameter space refinement) are provided below with reference to FIG. 8.

FIG. 8 illustrates steps of a method 800 for refining admissible parameter space of the simulation model 111 based on the failure region, according to an example embodiment. FIG. 8 is described below in conjunction with FIG. 7. In an example embodiment, the method 800 may be executed during the training of the failure classifier module 107b at step 703 of the calibration method 700. In another embodiment, the method 800 may be executed using a pre-trained failure classifier module 107b.

At step 801, the simulation model 111 (also referred to as “model”) corresponding to the industrial system 103 is obtained by the calibration system 101. In an example embodiment, the simulation model 111 may be predefined in the calibration system 101.

At step 803, failure of the simulation model 111 is defined. For example, the failure may be defined as unsuccessful simulation, forward in time, of the model 111 to a desired final time. In another embodiment, the failures may be defined as any simulation of the model 111 that takes more than a predetermined time interval for simulating the model 111. For example, if the model 111 is configured to complete the simulation within the predetermined time interval (for example 5 seconds), the failure is defined as any simulation that takes more than the predetermined time interval i.e., 5 seconds.

At step 805, the failure region is estimated using the FR-BO module 107, where the failure region is determined based on the definition of the failure defined at step 803. To that end, the FR-BO module 107 uses the failure classifier module 107b, where the failure classifier module 107b is trained to estimate the failure region based on values of the different combination of the parameters.

At step 807, the model 111 is refined i.e., the admissible parameter space of parameters in the different combinations of the parameters is refined based on the failure region estimated at step 805.

At step 809, it is determined whether the failure region of the refined model satisfies predetermined conditions for the failure region. The predetermined condition for the failure region may define size of the failure region within the admissible parameter space corresponding to the different combinations of the parameters. When it is determined that the failure region of the refined model does not satisfy the predetermined condition, for example, the failure region is not equal to or less than the predefined size of the failure region, the method 800 loops back, iteratively, to the step 807, where the refined model is refined again and the failure region is re-estimated.

On the other hand, when it is determined that the failure region of the refined model satisfies the predetermined condition, the method 800 proceeds to step 811, where the finalized refined model is used for calibrating the parameters of the model 111.

EMBODIMENTS

The above description provides exemplary embodiments only, and is not intended to limit the scope, applicability, or configuration of the disclosure. Rather, the above description of the exemplary embodiments will provide those skilled in the art with an enabling description for implementing one or more exemplary embodiments. Contemplated are various changes that may be made in the function and arrangement of elements without departing from the spirit and scope of the subject matter disclosed as set forth in the appended claims.

Specific details are given in the above description to provide a thorough understanding of the embodiments. However, understood by one of ordinary skill in the art can be that the embodiments may be practiced without these specific details. For example, systems, processes, and other elements in the subject matter disclosed may be shown as components in block diagram form in order not to obscure the embodiments in unnecessary detail. In other instances, well-known processes, structures, and techniques may be shown without unnecessary detail in order to avoid obscuring the embodiments. Further, like reference numbers and designations in the various drawings indicate like elements.

Also, individual embodiments may be described as a process which is depicted as a flowchart, a flow diagram, a data flow diagram, a structure diagram, or a block diagram. Although a flowchart may describe the operations as a sequential process, many of the operations can be performed in parallel or concurrently. In addition, the order of the operations may be re-arranged. A process may be terminated when its operations are completed but may have additional steps not discussed or included in a figure. Furthermore, not all operations in any particularly described process may occur in all embodiments. A process may correspond to a method, a function, a procedure, a subroutine, a subprogram, etc. When a process corresponds to a function, the function's termination can correspond to a return of the function to the calling function or the main function.

Furthermore, embodiments of the subject matter disclosed may be implemented, at least in part, either manually or automatically. Manual or automatic implementations may be executed, or at least assisted, through the use of machines, hardware, software, firmware, middleware, microcode, hardware description languages, or any combination thereof. When implemented in software, firmware, middleware or microcode, the program code or code segments to perform the necessary tasks may be stored in a machine readable medium. A processor(s) may perform the necessary tasks.

Various methods or processes outlined herein may be coded as software that is executable on one or more processors that employ any one of a variety of operating systems or platforms. Additionally, such software may be written using any of a number of suitable programming languages and/or programming or scripting tools, and also may be compiled as executable machine language code or intermediate code that is executed on a framework or virtual machine. Typically, the functionality of the program modules may be combined or distributed as desired in various embodiments.

Embodiments of the present disclosure may be embodied as a method, of which an example has been provided. The acts performed as part of the method may be ordered in any suitable way. Accordingly, embodiments may be constructed in which acts are performed in an order different than illustrated, which may include performing some acts concurrently, even though shown as sequential acts in illustrative embodiments. Although the present disclosure has been described with reference to certain preferred embodiments, it is to be understood that various other adaptations and modifications can be made within the spirit and scope of the present disclosure. Therefore, it is the aspect of the append claims to cover all such variations and modifications as come within the true spirit and scope of the present disclosure.

Claims

1. A calibration system for calibrating a model of dynamics of an industrial system that explains changes of a state of the industrial system in response to controlling the industrial system according to a sequence of control inputs to fit operational data indicative of measurements of an operation of the industrial system including values of the control inputs and corresponding values of the state of the industrial system, the calibration system comprising: at least one processor; and a memory having instructions stored thereon that, when executed by the at least one processor, cause the calibration system to: simulate an operation of the industrial system multiple times, each simulation cycle includes execution of the model of the industrial system having different combinations of values of parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation;train a failure classifier defining a likelihood of failure of the operation of the industrial system for the admissible range of values of parameters of the model using training data including the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation;train a parameter-to-cost regressor for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model of the industrial system and their corresponding calibration errors, wherein the parameter-to-cost regressor is trained iteratively until a termination condition is met, and wherein to perform an iteration the at least one processor is configured to: execute an acquisition function of a first two order moments of the calibration errors configured to identify a combination of parameters having maximum likelihood of minimizing the calibration errors according to the probabilistic parameter-to-cost mapping;update the identified combination of parameters with a failure-robust acquisition function adjusting the identified parameters according to the failure classifier to reduce the likelihood of the failure of the operation of the industrial system with the model having the updated parameters;simulate the operation of the industrial system controlled by the control inputs in the operational data with the model having the updated parameters to produce simulated states of the industrial system; andupdate the probabilistic parameter-to-cost mapping based on the updated parameters and calibration errors between the simulated states of the industrial system and corresponding states of the industrial system in the operational data; andcalibrate, when the termination condition is met, the model of the industrial system with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping according to the acquisition function.

2. The calibration system of claim 1, wherein the at least one processor is further configured to: estimate a failure region within the admissible range of values of parameters in the selected combination of parameters, wherein the failure region comprises at least one combination of values of the parameters that causes failure of the operation of the industrial system; andprovide probabilities of simulation failure over the entire admissible range of values of parameters in the selected combination of parameters based on the estimated failure region.

3. The calibration system of claim 1, wherein the training data further includes values of a cost function, wherein the cost function corresponds to the calibration error, and wherein the values of the cost function correspond to at least one of: real-valued scalers when the simulation of the model is successful or null values when the simulation of the model is unsuccessful.

4. The calibration system of claim 2, wherein, to obtain the failure region, the at least one processor is further configured to: sample within the admissible range of parameters to obtain the training data for learning the failure region, wherein the failure region is learnt using a scalable variational Gaussian process classifier (VGPC).

5. The calibration system of claim 4, wherein the parameter-to-cost regressor is trained iteratively using a Bayesian optimization, and wherein the Bayesian optimization includes a cost/reward function model using a Gaussian process (GP) for computing a probabilistic GP surrogate model for providing probabilistic parameter-to-cost mapping.

6. The calibration system of claim 4, wherein the at least one processor is further configured to: obtain predetermined conditions for the failure region; andestimate the failure region iteratively until the estimated failure region satisfies the predetermined condition for the failure region.

7. The calibration system of claim 5, wherein a prior function of the VGPC is initialized based on the training data and a user-specified kernel function.

8. The calibration system of claim 5, wherein the Bayesian optimization includes the acquisition function that exploits the probabilistic mapping provided by the probabilistic GP surrogate model to direct querying of consequent combination of the different parameters.

9. The calibration system of claim 5, wherein the Bayesian optimization uses a mean and a variance of the surrogate GP model to construct the acquisition function to select the combination of parameters having the maximum likelihood of minimizing the calibration error according to the probabilistic parameter-to-cost mapping.

10. The calibration system of claim 7, wherein the acquisition function includes an expected improvement (EI) acquisition function.

11. The calibration system of claim 7, wherein the user-specified kernel function comprises at least one of: a squared exponential kernel or a Matern Kernel.

12. The calibration system of claim 8, wherein the at least one processor is further configured to: obtain, using the scalable VPGC, probabilities of simulation failure over the entire admissible range of values of parameters in the selected combination of parameters, wherein the obtained probabilities of simulation failure are used by the acquisition function to form a failure-robust EI (FREI) acquisition function.

13. A calibration method for calibrating a model of dynamics of an industrial system that explains changes of a state of the industrial system in response to controlling the industrial system according to a sequence of control inputs to fit operational data indicative of measurements of an operation of the industrial system including values of the control inputs and corresponding values of the state of the industrial system, the calibration method comprising: simulating an operation of the industrial system multiple times, each simulation cycle includes execution of the model of the industrial system having different combinations of values of parameters selected within an admissible range of values of the parameters to estimate success or failure of the simulation;training a failure classifier defining a likelihood of failure of the operation of the industrial system for the admissible range of values of parameters of the model using training data including the selected combinations of the values of the parameters labeled with the estimated success or failure of the simulation;training a parameter-to-cost regressor for a probabilistic parameter-to-cost mapping between various combinations of different values of the parameters of the model of the industrial system and their corresponding calibration errors, wherein the parameter-to-cost regressor is trained iteratively until a termination condition is met, and wherein to perform an iteration the calibration method comprises: executing an acquisition function of a first two order moments of the calibration errors identifying a combination of parameters having the maximum likelihood of minimizing the calibration errors according to the probabilistic parameter-to-cost mapping;updating the identified combination of parameters with a failure-robust acquisition function adjusting the identified parameters according to the failure classifier to reduce the likelihood of the failure of the operation of the industrial system with the model having the updated parameters;simulating the operation of the industrial system controlled by the control inputs in the operational data with the model having the updated parameters to produce simulated states of the industrial system; andupdating the probabilistic parameter-to-cost mapping based on the updated parameters and calibration errors between the simulated states of the industrial system and corresponding states of the industrial system in the operational data; andcalibrating, when the termination condition is met, the model of the industrial system with an optimal combination of the parameters having the largest likelihood of minimizing the calibration errors at the probabilistic parameter-to-cost mapping according to the acquisition function.

14. The calibration method of claim 13, further comprising: estimating a failure region within the admissible range of values of parameters in the selected combination of parameters, wherein the failure region comprises at least one combination of values of the parameters that causes failure of the operation of the industrial system; andproviding probabilities of simulation failure over the entire admissible range of values of parameters in the selected combination of parameters based on the estimated failure region.

15. The calibration method of claim 14, wherein, for obtaining the failure regions the calibration, the calibration method further comprising: sampling within the admissible range of parameters to obtain the training data for learning the failure region, and wherein the failure region is learnt using a scalable variational Gaussian process classifier (VGPC).

16. The calibration method of claim 14, wherein the probabilistic parameter-to-cost regressor is trained iteratively using a Bayesian optimization, and wherein the Bayesian optimization includes a cost/reward function model using a Gaussian process (GP) for computing a probabilistic GP surrogate model for providing probabilistic parameter-to-cost mapping.

17. The calibration system of claim 14, wherein the method further comprises: obtaining predetermined conditions for the failure region; andestimating the failure region iteratively until the estimated failure region satisfies the predetermined condition for the failure region.