Learning Data-Driven Stable Corrections of Dynamical Systems. Application to the Simulation of the Top-Oil Temperature Evolution of a Power Transformer

1. Introduction

The two most usual approaches involved in simulation based engineering –SBE–, the one mainly based on physics, and the more recent one, the one based on the use of data manipulated by advanced machine learning techniques, both present their inherent limitations.

In general, the physics-based modeling framework produces responses that approximate quite well the real ones, as soon as an accurate enough model exists. The main difficulties when considering such a physics-based approach are: (i) the fidelity of the model itself, assumed calibrated; (ii) the impact that variability and uncertainty induces; and (iii) the computing time required for solving the complex and intricate mathematical models.

On the other hand, the data-driven framework is not fully satisfactory, because even if the data (assumed noise-free) represent well the reality, its extrapolation or interpolation in space and time from the collected data (in particular locations and times instants) entail usually a noticeable accuracy loss.

Of course, by increasing the number of collected data, one could expect approximating the real solution with more fidelity, however, data is not always simple to collect, not always possible to access, and in all cases collecting data is expensive (cost of sensors, cost of communication and analysis, ...). Equipping a very large industrial or civil infrastructure with millions of sensors to cover all its spatial dimension seems simply unreasonable.

Moreover, even when the solution is properly approximated, two difficulties persist: (i) the solution explainability, compulsory to certify solutions and decisions; and (ii) the domain of validity when extrapolating far from the domain where data was collected.

In view of the limitations of both existing procedures, a gateway consists of allying both to conciliate agility and fidelity. The hybrid paradigm seems a valuable and appealing option. It considers the reality expressible from the addition of two contributions: the existing knowledge (the state-of-the-art physics-based model or any other kind of knowledge-based model) and the part of the reality that the model ignores, the so-called ignorance (also called deviation, gap, discrepancy, or simply, ignorance).

1.1. The three main SBE methodologies revisited

To introduce and discuss different simulation based engineering frameworks, we consider a certain field $u\left(\mathbf{x}\right)$ describing the evolution of the variable u along the space defined by the coordinates $\mathbf{x}\in \Omega \subset {\mathbb{R}}^{\mathtt{D}}$.

We assume that it exists a certain knowledge on the addressed physics, described by the model $\mathcal{M}$, that in general consists of a system of algebraic equations or a system of differential equations complemented with the appropriate boundary and initial conditions ensuring the problem solvability. The solution provided by the just referred model, that as indicated, represents the existing knowledge on the problem at hand, is represented by $u^M\left(\mathbf{x}\right)$.

Due to the partial knowledge on the addressed physical phenomenon, the calculated solution $u^M\left(\mathbf{x}\right)$ is expected to differ from the reference one $u\left(\mathbf{x}\right)$, which is intended to be represented with maximum fidelity.

Thus, we define the residual $R^M\left(\mathbf{x}\right)$ according to

$$R^M\left(\mathbf{x}\right)=u\left(\mathbf{x}\right)-u^M\left(\mathbf{x}\right),$$

(1)

where the error can be computed from its norm, ${\mathcal{E}}^M=\parallel R^M\left(\mathbf{x}\right)\parallel$.

For reducing that error different possibilities exist:

• Physics-based model improvement. This approach consists of refining the modelling, by enriching the model itself, $\mathcal{M}\to \widehat{\mathcal{M}}$, such that its solution $u^{\widehat{M}}$ exhibits a smaller error, i.e. ${\mathcal{E}}^{\widehat{M}}\le {\mathcal{E}}^M$.
• Fully data-driven description. The data-driven route consists of making a large sampling of the space $\Omega$, ${\mathbf{x}}_1,...,{\mathbf{x}}_{\mathtt{P}}$, with $\mathtt{P}$ large enough, and with the points location ${\mathbf{x}}_i$, $i=1,...,\mathtt{P}$, maximizing the domain $\Omega$ coverage. These points are grouped into the set $\mathcal{S}$.
  The coverage is defined by the convex hull $\omega \left(\mathcal{S}\right)$ of the set $\mathcal{S}=\{{\mathbf{x}}_1,...,{\mathbf{x}}_{\mathtt{P}}\}$, ensuring interpolation for $\mathbf{x}\in \omega \left(\mathcal{S}\right)$, limiting the risky extrapolation to the region of $\Omega$ outside the convex hull $\omega \left(\mathcal{S}\right)$.
  Factorial samplings try to maximize the coverage, however factorial samplings, or the ones based on the use of Gauss-Lobatto quadratures, related to approximations making use of orthogonal polynomials [3], fail when the dimensionality of the space $\mathtt{D}$ increases.
  When $\mathtt{D}\gg 1$, sparse sampling is preferred, LHP (Latin hyper cube) for instance. Samplings based on gaussian processes –GP– aim at distributing the points in locations where the uncertainty is maximum (with respect to the predictions inferred from the previously collected data).
  Finally, the so-called active learning techniques drive the sampling aiming at maximizing the representation of a certain goal oriented quantity of interest [32].
  In what follows, we assume a generic sampling to access to the reference solution $u\left({\mathbf{x}}_i\right)$, $i=1,...,\mathtt{P}$, assuming a perfect measurability.
  Now, to infer the solution at $\mathbf{x}\in \omega \left(\mathcal{S}\right)$, it suffices constructing an interpolation or approximation, more generally, an adequate regression $u^{\mathcal{S}}\left(\mathbf{x}\right)$:

  $$u^{\mathcal{S}}\left(\mathbf{x}\right)\equiv u(\mathbf{x};u_1,...,u_{\mathtt{P}}),$$

  (2)

  with $u_1\equiv u\left({\mathbf{x}}_1\right),...,u_{\mathtt{P}}\equiv u\left({\mathbf{x}}_{\mathtt{P}}\right)$.
  Different possibilities exist, among them: regularized polynomial regressions Please rearrange all the references to appear in numerical order. [29], neural networks –NN– [16,31], support vector regression –SVR– [9], decision trees and their random forest counterparts [4,23], ... to name few.
  The trickiest issue concerns the error evaluation, that is quantified from a part of the data kept outside the training-set, the so-called test-set, used to quantify the performances of the trained regression.
  The main drawbacks of such a general procedure, particularly exacerbated in the multi-dimensional case ($\mathtt{D}\gg 1$), are the following:

  –

  Ability to explain the regression $u^{\mathcal{S}}\left(\mathbf{x}\right)$.

  –

  The size of the data-set ($\mathtt{P}$), scaling with the problem dimensionality $\mathtt{D}$.

  –

  The optimal sampling to cover $\Omega$ while guaranteeing the accuracy of $u^{\mathcal{S}}\left(\mathbf{x}\right)$, or the one of the goal oriented quantities of interest.

• Hybrid approach. It proceeds by embracing the physics-based and data-driven approaches, that as described in the next section, will improve the physics-based accuracy (while profiting of the physics-based explanatory capabilities) from the use of a data-driven enrichment, that at its turn, and under certain conditions, needs less amount of data than the fully data-driven approach just discussed [6].

2. Illustrating the hybrid approach

2.1. A simple linear regression reasoning

In the hybrid approach we consider first the contribution of a model expected representing reasonably the physics at hand, and noted by $u\left(\mathbf{x}\right)$, and that was noted previously as $u^M\left(\mathbf{x}\right)$, with $\mathbf{x}\in \Omega \subset {\mathbb{R}}^{\mathtt{D}}$.

As discussed before, it entails a residual $R^M\left(\mathbf{x}\right)$ and the associated error, ${\mathcal{E}}^M$.

Imagine for a while the simplest (and cheapest) regression of that residual, a linear regression involving $\mathtt{P}$ data, associated with the locations ${\mathbf{x}}_i\in \Omega ,i=1,...,\mathtt{P}$.

Thus, the linear regression involving a linear approximation of the residual, noted by $\Delta u^{\mathtt{L}}\left(\mathbf{x}\right)$ reads:

$$\Delta u^{\mathtt{L}}(\mathbf{x};a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta })=a_0^{\Delta }+a_1^{\Delta }{x}_1+...+a_{\mathtt{D}}^{\Delta }{x}_{\mathtt{D}},$$

(3)

that involve the unknown coefficients $a_0^{\Delta },a_1^{\Delta },...,a_{\mathtt{D}}^{\Delta }$.

By using the L2-norm, the linear regression results from the least-square minimization problem:

$$\{a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta }\}={\mathtt{argmin}}_{a_0^{*},...,a_{\mathtt{D}}^{*}}\left\{\sum _{i=1}^{\mathtt{P}}{\left(\Delta u^{\mathtt{L}}({\mathbf{x}}_i;a_0^{*},...,a_{\mathtt{D}}^{*})-R^M\left({\mathbf{x}}_i\right)\right)}^2\right\},$$

(4)

that defines an interpolation in the case $\mathtt{P}=\mathtt{D}+1$ and an approximation when $\mathtt{P}\ge \mathtt{D}+1$.

If we assume for a while that $u\left(\mathbf{x}\right)$ is fully available, and then $R^M\left(\mathbf{x}\right)$ too, the previous expression can be rewritten in a continuous form

$$\{a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta }\}={\mathtt{argmin}}_{a_0^{*},...,a_{\mathtt{D}}^{*}}{\parallel \Delta u^{\mathtt{L}}(\mathbf{x};a_0^{*},...,a_{\mathtt{D}}^{*})-R^M\left(\mathbf{x}\right)\parallel }_2.$$

(5)

The linear regression of the reference solution $u\left(\mathbf{x}\right)$, noted by $u^{\mathtt{L}}(\mathbf{x};a_0^u,...,a_{\mathtt{D}}^u)$, reads at its turn

$$u^{\mathtt{L}}(\mathbf{x};a_0^u,...,a_{\mathtt{D}}^u)=a_0^u+a_1^u{x}_1+...+a_{\mathtt{D}}^u{x}_{\mathtt{D}},$$

(6)

where the coefficients $a_0^u,...,u_{\mathtt{D}}^u$ results again from the least-square minimization problem

$$\{a_0^u,...,a_{\mathtt{D}}^u\}={\mathtt{argmin}}_{a_0^{*},...,a_{\mathtt{D}}^{*}}\left\{\sum _{i=1}^{\mathtt{P}}{\left(u^{\mathtt{L}}({\mathbf{x}}_i;a_0^{*},...,a_{\mathtt{D}}^{*})-u\left({\mathbf{x}}_i\right)\right)}^2\right\},$$

(7)

whose continuous expression reads

$$\{a_0^u,...,a_{\mathtt{D}}^u\}={\mathtt{argmin}}_{a_0^{*},...,a_{\mathtt{D}}^{*}}{\parallel u^{\mathtt{L}}(\mathbf{x};a_0^{*},...,a_{\mathtt{D}}^{*})-u\left(\mathbf{x}\right)\parallel }_2.$$

(8)

Eq. (5) implies

$$\parallel \Delta u^{\mathtt{L}}(\mathbf{x};a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta })-R^M\left(\mathbf{x}\right){\parallel }_2\le {\parallel R^M\left(\mathbf{x}\right)\parallel }_2.$$

(9)

Thus, as soon as the error related to the residual becomes smaller than the one related to the linear approximation of the reference solution, i.e. if

$${\mathcal{E}}^M=\parallel R^M{\left(\mathbf{x}\right)\parallel }_2\le {\parallel u^{\mathtt{L}}(\mathbf{x},a_0^u,...,a_{\mathtt{D}}^u)-u\left(\mathbf{x}\right)\parallel }_2,$$

(10)

then

$$\parallel \Delta u^{\mathtt{L}}(\mathbf{x};a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta })-R^M\left(\mathbf{x}\right){\parallel }_2\le \parallel R^M\left(\mathbf{x}\right){\parallel }_2\le {\parallel u^{\mathtt{L}}(\mathbf{x},a_0^u,...,a_{\mathtt{D}}^u)-u\left(\mathbf{x}\right)\parallel }_2,$$

(11)

that proves the higher accuracy of the hybrid approximation as soon as the solution provided by the model $\mathcal{M}$ represents a better approximation to the reference solution $u\left(\mathbf{x}\right)$ than a linear regression could attain.

2.2. General remarks

The analysis just addressed considers a simple linear regression involving a linear approximation, however it remains valid when considering regressions involving richer nonlinear approximations.

The first part of Eq. (11), i.e.

$$\parallel \Delta u^{\mathtt{L}}(\mathbf{x};a_0^{\Delta },...,a_{\mathtt{D}}^{\Delta })-R^M\left(\mathbf{x}\right){\parallel }_2\le {\parallel R^M\left(\mathbf{x}\right)\parallel }_2,$$

affirms that the simulated model solution is enriched in a L2-norm sense, and under the constraint of having enough amount of data to evaluate the considered norms. It is important to note, that even if the simulated model solution is enriched in the L2-norm sense, locally, the accuracy of the enriched solution could be degraded.

The second part pf Eq. (11), i.e.

$$\parallel R^M{\left(\mathbf{x}\right)\parallel }_2\le {\parallel u^{\mathtt{L}}(\mathbf{x},a_0^u,...,a_{\mathtt{D}}^u)-u\left(\mathbf{x}\right)\parallel }_2,$$

works under some conditions involving the data and the model, i.e. on $u\left(x\right)$ and $u^M\left(x\right)$, to be checked, and again needs having an amount of data enabling the calculation of the L2-norms.

2.3. On the domain of application of the hybrid modeling approach

The hybrid approach applies in different settings. It can be viewed as a sort of transfer learning, where rich behaviors are approached from others close enough and well stablished.

This hybrid approach seem particularly appealing in different situations, as the ones reported below:

• When the model captures most of the solution features, the correction must describe a discrepancy that exhibits smaller nonlinearities, as was the case treated in [27,28], where the same amount of data performed better within the hybrid than within the fully data-driven framework.
• Sometimes, the physics-based model operates very accurately in a part of the domain, whereas strong nonlinearities localize in a small region that can, in that case, be captured by a data-driven learned model, as considered in [26] for addressing the inelastic behavior of spot-welds.
• When considering the constitutive modeling of materials, the augmented rationale (or hybrid paradigm) expresses the real behavior from a first order behavior (calibrated from the available data) complemented by an enrichment (or correction) filling the gap between the collected data and the predictions obtained from the assumed model [15].
• When addressing plates (or shells) which noticeable 3D behaviors (deviating from the usual shell theory) a valuable solution consists of using an enriched kinematics consisting of two contributions: the first order one (usual shell kinematics) enriched with a second order contribution [25], that can be learned from data.
• The hybrid modeling can also transfer the existing knowledge slightly outside its domain of applicability with small amount of collected data, as performed in [24] for correcting state-of-the-art structural beam models.
• Sometimes the discrepancy concerns an imperfect alignment of the solution between the prediction and the measures. That discrepancy seems very high when evaluating it at each location, however, a small transport allows aligning both solutions. Optimal transport is very suitable in these situations, where the hybrid model consists of the usual nominal model enriched from a parametric correction formulated in an optimal transport setting as described in [34,35].
• In [20] a correction of a Mises yield function was performed from the deviation between the results predicted by using it, and the measures obtained at the structure level.
• Finally, when addressing processing and performances, the hybrid formulation, of different granularities, can exhibit advantages (amount of data, ability to explain, knowledge transfer, ...), at the heart of the digital twin developments [1,5,6,13,14,21,36].

3. Methods

In what follows, a system characterized by the state x that evolves in time, i.e. $x\left(t\right)$, is considered, and the stability of its numerical integration discussed [11]. When the state is multi-valued it will be noted by $\mathbf{x}$.

3.1. On the integration stability

We consider for a while a simple linear dynamical system expressed by

$$\frac{dx}{dt}=ax,a\in \mathbb{R}.$$

(12)

The integration reads

$$\frac{dx}{x}=adt\to \mathtt{ln}\left(x\right)=C+at\to x\left(t\right)=\widehat{C}{e}^{at},$$

(13)

with $\widehat{C}=e^C$ and $\widehat{C}=x(t=0)$.

It can be noticed that the existence of a bounded solution requires that $a\le 0$, because if $a>0$, $x(t\to \infty )=\infty$. Thus, the stability condition writes $a\le 0$.

Now, we consider the discrete case, in which the first order derivative is discretized by a first order finite difference

$$\frac{dx}{dt}=\frac{x_n-x_{n-1}}{\Delta t},$$

(14)

where $x_n=x(t=t_n)$, with $t_n=n\Delta t$.

In that case, the discrete time evolution reads

$$x_n=a{x}_{n-1}\Delta t+x_{n-1}=(a\Delta t+1)x_{n-1}=b{x}_{n-1},$$

(15)

with $b=(a\Delta t+1)$.

Now because $x_{n-1}=b{x}_{n-2}$, and so on, it finally results

$$x_n=b^n{x}_0,$$

(16)

from which it can be concluded that bounded solutions are subjected to the constraint $b\le 1$.

Similar results can be obtained in the case of multi-valued states. In what follows we consider the differential system

$$\frac{d\mathbf{x}}{dt}=\mathbf{A}\mathbf{x}.$$

(17)

For the sake of simplicity we assume that matrix $\mathbf{P}$ diagonalizes $\mathbf{A}$, with ${\mathbf{P}}^T\mathbf{P}=\mathbf{I}$ ($\mathbf{I}$ is the unit matrix). Thus, if we define $\mathbf{x}=\mathbf{P}\mathbf{y}$, it results

$$\mathbf{P}\frac{d\mathbf{y}}{dt}=\mathbf{A}\mathbf{P}\mathbf{y},$$

(18)

that premultiplying by ${\mathbf{P}}^T$ results

$${\mathbf{P}}^T\mathbf{P}\frac{d\mathbf{y}}{dt}={\mathbf{P}}^T\mathbf{A}\mathbf{P}\mathbf{y},$$

(19)

that taking into account that ${\mathbf{P}}^T\mathbf{P}=\mathbf{I}$, can be rewritten as

$$\frac{d\mathbf{y}}{dt}=\mathbf{D}\mathbf{y},$$

(20)

with $\mathbf{D}$ diagonal.

The fact of being $\mathbf{D}$ diagonal, allows decoupling the solution of Eq. (20). Thus, we have

$$\frac{d{y}_i}{dt}=D_i{y}_i,\forall i,$$

(21)

with $y_i$ the i-component of $\mathbf{y}$ and $D_i$ the $(i,i)$-diagonal component of $\mathbf{D}$.

Thus, using the previous results, the stability condition reads

$${\mathtt{max}}_i{D}_i\le 0.$$

(22)

By discretizing Eq. (20), it results

$${\mathbf{y}}_n=(\Delta t\mathbf{D}+\mathbf{I}){\mathbf{y}}_{n-1},$$

(23)

or notting $\widehat{\mathbf{D}}=\Delta t\mathbf{D}+\mathbf{I}$,

$${\mathbf{y}}_n=\widehat{\mathbf{D}}{\mathbf{y}}_{n-1},$$

(24)

that using the recurrence results in

$${\mathbf{y}}_n={\widehat{\mathbf{D}}}^n{\mathbf{y}}_0.$$

(25)

Thus, stability requires the spectral radius $\rho$ of $\widehat{\mathbf{D}}$ being lower than the unity, i.e. $\rho \left(\widehat{\mathbf{D}}\right)\le 1$.

These conditions should be satisfied by the dynamical learned models, as discussed later.

3.2. Learning integrators

When the dynamical system is know, different numerical integration schemes (explicit, implicit or semi-implicit) can be applied, with different convergence orders, for instance the Euler, or the more accurate Runge-Kutta, among many other choices. Thus, the discretization time-step $\Delta t$ must be chosen to ensure stability, and to guarantee convergence, that is, a numerical solution close enough to the reference solution of the problem.

In what follows, we revisit, and discuss, different machine learning techniques that enable learning dynamical systems and applying them for the integration of dynamical systems, as recurrent NN –rNN– and long short time memory –LSTM– techniques perform [10,37].

In this section $\mathbf{h}$ refers to the hidden state, whereas $\mathbf{l}$ and $\mathbf{o}$ refer, respectively, to the inputs (loading) and observable outputs.

3.2.1. Recurrent Neural Network

If ${\mathbf{W}}^{•}$ represents dense matrices associated to variables •, $\sigma (·)$ the activation function, ${\mathbf{h}}_n$ the internal state at the n-time step, $t_n$, and ${\mathbf{l}}_n$ and ${\mathbf{o}}_n$ the associated input and output at the present time $t_n$, then the rNN proceeds from:

$$\left\{\begin{array}{c}{\mathbf{h}}_n={\sigma }^h({\mathbf{W}}^h{\mathbf{h}}_{n-1}+{\mathbf{W}}^l{\mathbf{l}}_n)\hfill \\ {\mathbf{o}}_n={\sigma }^y\left({\mathbf{W}}^o{\mathbf{h}}_n\right),\hfill \end{array}\right.,$$

(26)

whose architecture is depicted in Figure 1.

To perform the time evolution it suffices to use the rNN and re-inject the rNN output ${\mathbf{h}}_n$, as the input at the next time step.

A temporal sequence of inputs and outputs, assumed available, allows calculating the different weights composing the different $\mathbf{W}$ matrices.

To highligth the connexion between the rNN and usual dynamical system integrators, we consider for a while the finite element semi-discretized form of a linear parabolic differential equation

$$\mathbf{M}\frac{d\mathbf{h}}{dt}=\mathbf{K}\mathbf{h}+\mathbf{l},$$

(27)

where $\mathbf{M}$ and $\mathbf{K}$ are two matrices that results from the spatial discretization. Here, $\mathbf{h}$ refers to the state and $\mathbf{l}$ to the system loading.

Then, the explicit discretization of Eq. (27) reads

$$\mathbf{M}{\mathbf{h}}_n=\Delta t\mathbf{K}{\mathbf{h}}_{n-1}+\Delta t{\mathbf{l}}_n+\mathbf{M}{\mathbf{h}}_{n-1}=(\Delta t\mathbf{K}+\mathbf{M}){\mathbf{h}}_{n-1}+\Delta t{\mathbf{l}}_n,$$

(28)

from which the state updating results

$${\mathbf{h}}_n={\mathbf{M}}^{-1}(\Delta t\mathbf{K}+\mathbf{M}){\mathbf{h}}_{n-1}+\Delta t{\mathbf{M}}^{-1}{\mathbf{l}}_n,$$

(29)

that can be rewritten as

$${\mathbf{h}}_n={\mathbf{W}}^h{\mathbf{h}}_{n-1}+{\mathbf{W}}^l{\mathbf{l}}_n,$$

(30)

that correspond to the particularization of Eq. (26) to the linear case. Thus, the intimate connexion between rNN and usual discretization techniques becomes explicit. The last operates on a known model whereas the former learns the model itself (the different ${\mathbf{W}}^{•}$ matrices).

Remark 1. In the just discussed linear model, when the whole state is observed, i.e. ${\mathbf{o}}_n={\mathbf{h}}_n$, ${\mathbf{W}}^o=\mathbf{I}$.

3.2.2. On the model memory

When considering a first order dynamical system, as just discussed, the solution at the present time can be computed from the only knowledge of the previous solution (at the previous time step) and the present loading (also known as action or input). Thus, one could imagine that for learning first order dynamical systems a short memory, as for instance the one of the just described rNN, suffices.

However, sometimes, larger memory are needed even when addressing first order models as described in this section.

For the sake of simplicity we consider the system state given by $\mathbf{h}\left(t\right)={(h_1\left(t\right),h_2\left(t\right))}^T$, whose evolution is governed by a simple linear first order dynamical system

$$\left\{\begin{array}{c}{\dot{h}}_1\left(t\right)=K_{11}{h}_1\left(t\right)+K_{12}{h}_2\left(t\right)+l_1\left(t\right)\hfill \\ {\dot{h}}_2\left(t\right)=K_{21}{h}_1\left(t\right)+K_{22}{h}_2\left(t\right)+l_2\left(t\right)\hfill \end{array}\right.,$$

(31)

where, without loss of generality, it is assumed $K_{22}=0$ and $l_2\left(t\right)=0$, i.e.

$$\left\{\begin{array}{c}{\dot{h}}_1\left(t\right)=K_{11}{h}_1\left(t\right)+K_{12}{h}_2\left(t\right)+l_1\left(t\right)\hfill \\ {\dot{h}}_2\left(t\right)=K_{21}{h}_1\left(t\right)\hfill \end{array}\right..$$

(32)

In what follows, it is assumed that $h_1\left(t\right)$ and $x_1\left(t\right)$ are accessible, but that nothing, even the existence of $l_2\left(t\right)$, is known. In that circumstances the key question is: can we learn a model relating $l_1\left(t\right)$ and $h_1\left(t\right)$, both accessible and measurable, knowing that the state $h_1\left(t\right)$ depends on another one, $h_2\left(t\right)$, that evolves in time and remains inaccessible, and consequently unknown.

To facilitate the model manipulation, the Fourier transform applies, here noted by the symbol ${·}^{*}$, acting on the different time-dependent variables. Thus, Eq. (32) writes

$$\left\{\begin{array}{c}i\omega h_1^{*}=K_{11}{h}_1^{*}+K_{12}{h}_2^{*}+l_1^{*}\hfill \\ i\omega h_2^{*}=K_{21}{h}_1^{*}\hfill \end{array}\right..$$

(33)

From the second equation in (33) we obtain

$$h_2^{*}=K_{21}\frac{h_1^{*}}{i\omega }$$

(34)

that inserted in the first equation in (33) leads to

$$i\omega h_1^{*}=K_{11}{h}_1^{*}+K_{12}{K}_{21}\frac{h_1^{*}}{i\omega }+l_1^{*},$$

(35)

that coming back to the time domain results

$$\frac{d{h}_1\left(t\right)}{dt}=K_{11}{h}_1\left(t\right)+K_{12}{K}_{21}{\int }_0^t{h}_1\left(\tau \right)d\tau +l_1\left(t\right),$$

(36)

that proves that ignoring components of the state manifests in the measurable state history, here symbolized by the integral. In a certain way, this result can be interpreted as a consequence of the Taken delay embedding theorem.

Thus, memory is not only a consequence of the order of the subjacent dynamics, the memory in the learning procedure also depends on the ignored states. Thus, sometimes, when addressing poorly described physical systems or partially accessible, larger memory than the one that recurrent NN offers, seems compulsory.

3.2.3. Learners with larger memory

LSTM (Long Short Time Memory) [19] becomes an appealing technique for learning time-dependent problems while ensuring larger memory. LSTM combines a short and a long time memory units, with an evanescent memory for the long-time path, and a combination of long and short leads to the short memory response..

A particular architecture is depicted in Figure 2, that is considered for including not only memory, but also the previous inputs as well as forgetting the initial conditions of the long ${\mathbf{c}}_0$ and short ${\mathbf{h}}_0$ memory channels, to facilitate its use as a generic time integrator.

Recurrent NN (rNN) and LSTM learn the state evolution, and special care must be paid to ensure time integration stability.

Sometimes it seems preferably learning the forcing term of the dynamical system, as Residual Nets or Neural Differential Equations techniques perform. Next section revisits Residual Nets and discuss stability issues.

3.2.4. Residual Nets

Residual Neural Networks (ResNet) [2,17,18] aim at emulating the backward integration of a dynamical system. When we consider the dynamical system

$$\frac{d\mathbf{x}}{dt}=\mathbf{f}\left(\mathbf{x}\right),\mathbf{x}(t=0)={\mathbf{x}}_0,$$

(37)

the ResNet looks for the function $f\left(\mathbf{x}\right)$ that allows reproducing the state time evolution $\mathbf{x}\left(t\right)$, from its discrete knowledge, that is, from ${\mathbf{x}}_n=\mathbf{x}\left(t_n\right)$, with $t_n=n\Delta t$.

By considering the simplest time stepping, the time derivative can be approximated by $d\mathbf{x}/dt=({\mathbf{x}}_{n+1}-{\mathbf{x}}_n)/\Delta t$, and then, one can use the following updating

$${\mathbf{x}}_{n+1}={\mathbf{x}}_n+\Delta t\mathbf{f}\left({\mathbf{x}}_n\right),$$

(38)

and learn the forcing term $\mathbf{f}\left(\mathbf{x}\right)$ that induced the state evolution by using an appropriate regression.

For that purpose ResNet creates a NN trained from ${\mathbf{x}}_n$ and ${\mathbf{x}}_{n+1}-{\mathbf{x}}_n$, $n=1,...$. Then, as soon as the NN is trained, the output ${\mathbf{x}}_{n+1}$ is calculated by applying on ${\mathbf{x}}_n$ the NN just constructed adding the identity operator applied on the input.

Now, the dynamical system integration consists in applying in a recursive way the ResNet, with the output at time n becoming the input of the next time-step. The ResNet architecture and the integration procedure are illustrated in Figure 3.

4. Simple procedures for stabilizing ResNets-based integration

In what follows, the stabilization of the time integration performed with ResNet is discussed. First, a simple stabilization is accomplished by linearizing the dynamical system, while using of the results just discussed.

Then, inspired by that accomplishment, a more general stabilization, able to operate in nonlinear settings, is proposed.

Finally, the performances of both strategies are compared.

4.1. Learning stable linear dynamical systems

In what follows we consider a parametrized dynamical system,

$$\frac{dx}{dt}=f(x;\mathbf{z}),$$

(39)

where the state at time $t_n$, $x_n$, depends on the previous state $x_{n-1}$ and the loading $f(·)$ driving the state change. The last is assumed depending not only on the state, but also on a series of parameters (input features) here grouped in vector $\mathbf{z}$. Thus, the learning procedure aims at computing the regression $f(x;\mathbf{z})$.

However, for obtaining a stable integrator, certain constraints must be fulfilled. As these constraints have an easy form in the linear case, as previously discussed, here we proceed to linearize the forcing term as follows

$$f(x;\mathbf{z})\approx g\left(\mathbf{z}\right)x+h\left(\mathbf{z}\right).$$

(40)

Now, the stability is ensured as soon as $g\left(\mathbf{z}\right)\le 0$, constraint that can be easily enforced by employing a penalty term in the loss function of the NN associated with $g\left(\mathbf{z}\right)$.

The linearized ResNet resembles to the DMD (dynamic mode decomposition) [30], with probably an easier introduction of the parametric dimension.

4.2. Learning stable nonlinear dynamical systems

Inspired from the just discussed rationale, a possible route for enforcing stability without detriment of the nonlinear behavior consists of writing

$$f(x;\mathbf{z})\approx g(\mathbf{z},x)x+h\left(\mathbf{z}\right).$$

(41)

while enforcing in the construction of regression $g(x;\mathbf{z})$ the negativity constraint, that is, $g(x;\mathbf{z})\le 0$.

4.3. Numerical examples and discussion

To prove the expected performances we consider two cases, one linear, where both methodologies just described should work, and a nonlinear one, where the former is expected failing, whereas the later (the nonlinear counterpart) is expected performing correctly.

This section only addresses the stability concerns, without giving details on the neural architectures employed for learning functions $g\left(\mathbf{z}\right)$, $\mathbf{g}(x;\mathbf{z})$ and $h\left(\mathbf{z}\right)$, that will be described in detail in Section 5. Moreover, when the learnt functions are employed for integrating the dynamical system incrementally, by using a standard first order explicit time-marching, the computed solution will be noted with hat symbol, $\widehat{·}$.

4.3.1. Linear dynamical system

We consider data generated from the integration of the dynamical system

$$\frac{dx}{dt}=1-x,$$

(42)

integrated from the initial condition $x(t=0)$= 0.

The long-time solution $x=1$ is stable, because if $x<1$ the solution increases, whereas if $x>1$, it decreases.

The solution reads

$$x\left(t\right)=1-e^{-t},$$

(43)

that is used for generating the synthetic data that will be used then, for learning the dynamical system by using the two procedures previously described.

Figure compares the solution obtained by integrating numerically the models learnt by employing both, the linear and nor linear learning procedures. The problem being linear both procedures are expected performing identically, as Figure proves.

4.3.2. Nonlinear dynamical system

Now, a nonlinear dynamical system is considered, expressed by

$$\frac{dx}{dt}={(1-x)}^2,$$

(44)

integrated from the initial condition $x(t=0)$= 0.

Using the same rationale than before, we can conclude on the stability of the associated long-time solution $x=1$.

Figure 5 compares the solution obtained with the linear and nonlinear learning procedures. As expected, the nonlinear learning procedure ensures stability and represents accurately the nonlinear behavior, whereas the learnt linear model results stable but remains inaccurate.

4.4. Final remarks

Another potential neural architecture with improved robustness with respect to stability is the NeuralODE [7,8,12], that learns the model while considering many integration steps, instead of learning for one-step time increments as ResNet performs. Even is stability is significantly enhanced, the NeuralODE training efforts are higher than the ones required by the ResNet, mainly due to the fact that the back-propagation in the loss function minimization needs the solution of an adjoint problem, that, in some cases exhibits poor convergence.

5. Application to the evaluation of the top-oil temperature of an electric power transformer

This section addresses an applicative case of industrial relevance. Aging of transformers exhibits a correlation with the temperature of the oil all along their lives in service, and moreover, the oil temperature seems to be an appealing feature to anticipate faults, and consequently to be used in predictive maintenance.

Some simplified models exist to evaluate the oil temperature depending on the transformer delivered power and the ambient temperature. Standards propose different simplified models for that purpose, as the one considered later.

Because of the complexity of a transformer, large in size and embracing numerous physics and scales, an appealing modeling route consists of using a simplified model and then enrich it from the available collected data [22].

The just referred correction comes from the time-integration of a dynamical system, involving a parametric loading (delivered power and ambient temperature), learnt from the available data under the stability constraints.

To illustrate the construction of the enriched model for the top oil temperature prediction $\Theta$, we consider as input the ambient temperature $T^{amb}$, and the transformer load $K^{load}$, which are both measured every hour. Thus, the model parameters read $\mathbf{z}={(T^{amb};K^{load})}^T$.

The transformer oil temperature results from

$$\Theta \left(t\right)={\Theta }^{TO}\left(t\right)+d\left(t\right),$$

(45)

where $d\left(t\right)$ represents the deviation, and ${\Theta }^{TO}$ the temperature predicted by a state-of-the-art simplified model [33]:

$$\left\{\begin{array}{c}P(K^{load},{\Theta }^W,X^{tp})=C_{th}\frac{d\Delta {\Theta }^{TO}}{dt}+\frac{\Delta {\Theta }^{TO}}{R_{th}}\hfill \\ {\Theta }^{TO}=T^{amb}+\Delta {\Theta }^{TO}\hfill \\ {\Theta }^W={\Theta }^{TO}+\Delta {\Theta }^{OW}\hfill \end{array}\right.$$

(46)

where

• $X^{tp}$: The position of the tap-changer;
• $K^{load}$: the load factor, a ratio between the nominal load current and the actual current;
• P: represents the iron, copper and supplementary losses. The power that heats the oil is composed of the losses that do not depend on the transformer load (iron losses, supposed constant) and the losses that depend on the transformer load (copper and supplementary losses) that depend on the average windings temperature and the load factor, according to: $P=P_0+P_{load}$, with $P_{load}=K{{}^{load}}^2{P}_{loa{d}_r}$ and $P_{loa{d}_r}=k{P}_{Joul{e}_r}+\frac{P_{sup{p}_r}}{k}$ (k being a correction factor related to the material resistivity);
• $T^{amb}$: The ambient temperature;
• ${\Theta }^{TO}$: The simulated top oil temperature
• $\Delta {\Theta }^{TO}$: The temperature difference between the simulated top oil temperature and the ambient temperature;
• $R_{th}$ and $C_{th}$: thermal resistance and thermal capacitance of the equivalent transformer thermal circuit;
• ${\Theta }^W$: The average winding temperature;
• $\Delta {\Theta }^{OW}$: The difference between the average winding temperature and the simulated oil temperature ${\Theta }^{TO}$. It is supposed constant and found during the commissioning test (standards)

The physics-based model (46) is calibrated from the available experimental data, from which, parameters $C_{th}$ and $R_{th}$ are obtained. Two physics-based models are used, a linear one where $C_{th}$ and $R_{th}$ are assumed constant, and a nonlinear physics-based model where $C_{th}$ and $R_{th}$ are temperature dependent, that is, both coefficients depend on the top-oil temperature ${\Theta }^{TO}$.

The available experimental data $\Theta \left(t\right)$, the prediction from the calibrated physics-based non-linear model (46), ${\Theta }^{TO}\left(t\right)$, and the deviation between both them, $d\left(t\right)=\Theta \left(t\right)-{\Theta }^{TO}$, are all depicted in Figure 6.

The model correction (the deviation model, $d\left(t\right)$) is obtained also using two different approaches: the linearized ResNet:

$$\frac{\partial d}{\partial t}=g\left(\tilde{\mathbf{z}}\right)d+h\left(\tilde{\mathbf{z}}\right),$$

(47)

and the nonlinear ResNet counterpart:

$$\frac{\partial d}{\partial t}=g(\tilde{\mathbf{z}},d)d+h\left(\tilde{\mathbf{z}}\right),$$

(48)

where $\tilde{\mathbf{z}}$ represents the augmented features set, that contains the physical features $\mathbf{z}$ augmented with the model prediction ${\Theta }^{TO}$.

Functions $g\left(\tilde{\mathbf{z}}\right)$ and $h\left(\tilde{\mathbf{z}}\right)$ are described by using two LSTM architectures (described in Table 1 and Table 2 respectively for the linearized ResNet) that consider the extended features involved in $\tilde{\mathbf{z}}$ at the present time, as well as at the previous 4 time-steps.

Thus, the linearized correction dynamical model reads

$$\left\{\begin{array}{c}{d}_n=\Delta tg\left(\tilde{\mathbf{z}}\right)d_{n-1}+\Delta th\left(\tilde{\mathbf{z}}\right)+d_{n-1}\hfill \\ \tilde{\mathbf{z}}=\left(\begin{array}{ccc}{T}_{n-4}^{amb}& K_{n-4}^{load}& {\Theta }_{n-4}^{TO}\\ T_{n-3}^{amb}& K_{n-3}^{load}& {\Theta }_{n-3}^{TO}\\ & ⋮& \\ T_n^{amb}& K_n^{load}& {\Theta }_n^{TO}\end{array}\right)\hfill \end{array}\right..$$

(49)

The neural network architecture considered for describing functions $g\left(\tilde{\mathbf{z}}\right)$ and $h\left(\tilde{\mathbf{z}}\right)$ are both based on the use of LSTM layers combined with a deep dense neural network layer, as described in Table 1 and Table 2 for the linearized ResNet. They were built-up by using Tensorflow Keras libraries. The inputs involved in Eq. (49) are shown in Figure 7.

The training was performed on the initial $80\%$ of the available measures, while the testing was performed on the remaining $20\%$.

The prediction performed from the corrected (enriched) model is depicted in Figure 8, on the testing interval. It is important to note, that the learnt models $g\left(\tilde{\mathbf{z}}\right)$ and $h\left(\tilde{\mathbf{z}}\right)$, were used to integrate the dynamical problem that governs the time evolution of the deviation, Eq. (49), that is, the computed deviation at time $t_i$ was computed from the one at time $t_{i-1}$, and then it served to compute the one at the next time step $t_{i+1}$, and so on.

To distinguish between the known deviation $d\left(t\right)$ and the one computed from the integrator based on the learnt dynamics $g\left(\tilde{\mathbf{z}}\right)$ and $h\left(\tilde{\mathbf{z}}\right)$, the later will be noter by the hat symbol, i.e.

$${\widehat{d}}_{n+1}=\Delta tg\left(\tilde{\mathbf{z}}\right){\widehat{d}}_n+\Delta th\left(\tilde{\mathbf{z}}\right)+{\widehat{d}}_n$$

(50)

From the computed correction, the model was enriched, exhibiting excellent accuracy and the stability performances.

The fact that the linearized dynamical system operating on the solution correction exhibits good performances, proves that most of the problem nonlinearities were captured by the first order simplified model. When it comes to the nonlinear version of the correction effort, the dynamic version of the model is written in a similar manner as shown in equation (49), with the dynamical integration form used:

$${\widehat{d}}_{n+1}=\Delta tg(\tilde{\mathbf{z}},{\widehat{d}}_n){\widehat{d}}_n+\Delta th\left(\tilde{\mathbf{z}}\right)+{\widehat{d}}_n$$

(51)

To compare the performances of the hybrid model (data-driven enrichment of the simplified physics-based model), the model governing the experimental data evolution was learnt using the same rationale but now applied on the data, according to:

$$\left\{\begin{array}{c}{\Theta }_{n+1}=\Delta t{g}^{\theta }\left(\mathbf{z}\right){\Theta }_n+\Delta t{h}^{\theta }\left(\mathbf{z}\right)+{\Theta }_n\hfill \\ \mathbf{z}=\left(\begin{array}{cc}{T}_{n-4}^{amb}& K_{n-4}^{load}\\ T_{n-3}^{amb}& K_{n-3}^{load}\\ & ⋮& \\ T_n^{amb}& K_n^{load}\end{array}\right)\hfill \end{array}\right.,$$

(52)

where $g^{\theta }\left(\mathbf{z}\right)$ and $h^{\theta }\left(\mathbf{z}\right)$ refer to the parametric functions related to the measured oil temperature, both them depending exclusively on the input features $\mathbf{z}$.

Again, the training was performed using the same model presented in Table 1 and Table 2, with the same $80\%$ and $20\%$ training and test data-sets.

Figure 9 depicts the results obtained from the integration, where again the hat symbol refers the integrated temperature.

$${\widehat{\Theta }}_{n+1}=\Delta t{g}^{\theta }\left(\mathbf{z}\right){\widehat{\Theta }}_n+\Delta t{h}^{\theta }\left(\mathbf{z}\right)+{\widehat{\Theta }}_n.$$

(53)

Figure 9 depicts the data-driven model integration, proving the stability performances. However, when comparing the residual error (with respect to the experimental data) of the fully data-driven prediction and the one related to the physics-based enriched solution obtained from the data-driven model of the deviation, both compared in Figure 10, the hybrid modeling framework performs better, ensuring higher accuracy.

Table 3, compares the mean value of the errors associated with the different tested models. From Table 3, one can infer that:

• The hybrid approach improves the physics-based model performances.
• Enriching a richer nonlinear physics-based model outperform with respect to enriching the linear counterpart of the simplified physics-based model.
• When the considered physics-based models are too far from the reference solution (experimental data) the data-driven model can outperform the hybrid modeling.

Table 3. Comparing different built models: mean error (in °C) on the testing set. The nonlinear physics-based model mean error was 3.91°C when used alone, over the same testing set, and 3.25°C for the linear physics based model.

ResNet	Fully	HT from a linear	HT from a nonlinear
	data-driven	physical model	physical model
Linear stabilized	2.173	3.143	1.620
Nonlinear stabilized	1.716	1.516	1.439

Table 3. Comparing different built models: mean error (in °C) on the testing set. The nonlinear physics-based model mean error was 3.91°C when used alone, over the same testing set, and 3.25°C for the linear physics based model.

ResNet	Fully	HT from a linear	HT from a nonlinear
	data-driven	physical model	physical model
Linear stabilized	2.173	3.143	1.620
Nonlinear stabilized	1.716	1.516	1.439

To show the effect of the stabilization, a ResNet is trained without enforcing the stability constraints previously proposed, by using the formulation:

$$\frac{\partial d}{\partial t}=g(\tilde{\mathbf{z}},d)+h\left(\tilde{\mathbf{z}}\right),$$

(54)

with no conditions imposed during the calculation of g and h.

The discrete form

$${\widehat{d}}_{n+1}=\Delta tg(\tilde{\mathbf{z}},d_n)+\Delta th\left(\tilde{\mathbf{z}}\right)+d_n,$$

(55)

provides excellent predictions. It is important note that here the solution at time $t_{n+1}$, ${\widehat{d}}_{n+1}$, is computed from the exact deviation at time $t_n$, $d_n$.

However, a full integration where ${\widehat{d}}_{n+1}$ is computed from the previously computed ${\widehat{d}}_n$

$${\widehat{d}}_{n+1}=\Delta tg(\tilde{\mathbf{z}},{\widehat{d}}_n)+\Delta th\left(\tilde{\mathbf{z}}\right)+{\widehat{d}}_n,$$

(56)

produces extremely bad predictions, a direct consequence of the lack of stability.

Figure 11 proves the importance of using stable formulations when the learnt model is expected serving for performing integrations from an initial condition, as it is always the case in prognosis applications.

6. Conclusion

The present paper proposed a hybrid framework where the data-driven model serves to enrich a physics-based model considered as a first approximation of the addressed physics.

We proved that the hybrid framework enhances the predictions accuracy with respect to the physics-based model, however, the hybrid approach is superior to a fully data-driven model only under certain conditions.

The learning technique employed to model the time evolution of the deviation (or the one of the data) must ensure the integration stability. A stabilization has been proposed and its performance proved.

The application to a problem of practical relevance proved the excellent performances of the proposed methodology.