Physics-informed State-space Neural Networks for Transport Phenomena
This work introduces Physics-informed State-space neural network Models
(PSMs), a novel solution to achieving real-time optimization, flexibility, and
fault tolerance in autonomous systems, particularly in transport-dominated
systems such as chemical, biomedical, and power plants. Traditional data-driven
methods fall short due to a lack of physical constraints like mass
conservation; PSMs address this issue by training deep neural networks with
sensor data and physics-informing using components' Partial Differential
Equations (PDEs), resulting in a physics-constrained, end-to-end differentiable
forward dynamics model. Through two in silico experiments -- a heated channel
and a cooling system loop -- we demonstrate that PSMs offer a more accurate
approach than a purely data-driven model. In the former experiment, PSMs
demonstrated significantly lower average root-mean-square errors across test
datasets compared to a purely data-driven neural network, with reductions of 44
%, 48 %, and 94 % in predicting pressure, velocity, and temperature,
respectively.
  Beyond accuracy, PSMs demonstrate a compelling multitask capability, making
them highly versatile. In this work, we showcase two: supervisory control of a
nonlinear system through a sequentially updated state-space representation and
the proposal of a diagnostic algorithm using residuals from each of the PDEs.
The former demonstrates PSMs' ability to handle constant and time-dependent
constraints, while the latter illustrates their value in system diagnostics and
fault detection. We further posit that PSMs could serve as a foundation for
Digital Twins, constantly updated digital representations of physical systems.
Akshay J. Dave, Richard B. Vilim
# Physics-informed State-space Neural Networks for Transport Phenomena

###### Abstract

This work introduces Physics-informed State-space neural network Models (PSMs), a novel solution to achieving real-time optimization, flexibility, and fault tolerance in autonomous systems, particularly in transport-dominated systems such as chemical, biomedical, and power plants. Traditional data-driven methods fall short due to a lack of physical constraints like mass conservation; PSMs address this issue by training deep neural networks with sensor data and physics-informing using components' Partial Differential Equations (PDEs), resulting in a physics-constrained, end-to-end differentiable forward dynamics model. Through two in silico experiments -- a heated channel and a cooling system loop -- we demonstrate that PSMs offer a more accurate approach than purely data-driven models.

Beyond accuracy, there are several compelling use cases for PSMs. In this work, we showcase two: the creation of a nonlinear supervisory controller through a sequentially updated state-space representation and the proposal of a diagnostic algorithm using residuals from each of the PDEs. The former demonstrates the ability of PSMs to handle both constant and time-dependent constraints, while the latter illustrates their value in system diagnostics and fault detection. We further posit that PSMs could serve as a foundation for Digital Twins, constantly updated digital representations of physical systems.

State-space model Physics-informed Data-driven

Supervisory Control Diagnostics

## 1 Introduction

Autonomous systems are envisioned to achieve real-time optimization, be flexible in operation, and be fault-tolerant [1]. These capabilities necessitate a model-based approach. However, purely data-driven methods, which have proliferated in the literature, lack physical constraints like mass conservation, precluding their adoption by industry. In this work, we propose a method that _fuses_ physical measurements from sensors with physical knowledge of the conservation laws. The application of this method is specifically developed for systems dominated by transport phenomena.

Transport phenomena, which describe the movement of mass, momentum, and energy, are fundamental to many engineering fields as they offer a predictive understanding of systems' behavior. In chemical engineering, these phenomena facilitate the design and optimization of processes in reactors, separators, and heat exchangers, which are essential for chemical, fuel, and pharmaceutical production [2]. Biomedical engineering applies transport phenomena to comprehend nutrient, oxygen, and waste transportation in biological systems like the circulatory system, driving the development of new medical technologies such as artificial organs [3]. Power plant engineering utilizes transport phenomena in optimizing heat and mass transfer in boilers and heat exchangers for efficient and safe power generation and in emission control [4]. They are also vital in nuclear reactor design, where effective heat transfer from the
reactor core to the coolant is paramount [5]. Consequently, modeling transport phenomena is crucial for the design, optimization, and enhancement of engineering systems and processes.

A set of transport equations models transport phenomena. The solution of transport equations involves formulating functions of space and time for multiple fields (e.g., pressure, velocity, temperature), given a set of boundary and initial conditions. Classically, there have been two major approaches to solving transport equations: analytical and numerical. The former applies to situations where the problem setting or transport phenomena are simple enough to directly use mathematical methods (e.g., separation of variables) - this is generally not the case for systems of interest described above. The latter is used for systems with complex geometries, tightly coupled fields, time-dependent boundary/initial conditions, or parameters that vary with the field itself. Numerical methods include finite difference, finite volume, or finite element [6]. These methods discretize the domain into a grid, and values at fixed locations represent the unknown fields. The transport equations are then approximated using numerical algorithms (e.g., the forward Euler method, the Crank-Nicolson method) and solved iteratively until convergence is achieved.

The numerical methods are the current state-of-the-art approach to modeling transport phenomena for engineering systems. Their main advantage is versatility in addressing a range of complex problems with time-dependent boundary/initial conditions. However, they have several disadvantages. _Discretization error_: numerical methods discretize the domain and approximate the transport equations. This process introduces error, which can be reduced with finer grids - unfortunately increasing computational cost. _Convergence_: because numerical methods involve iterative algorithms, convergence to a solution is not guaranteed and can result in numerical instability. _Static assumptions_: generally, several parameters in a transport equation must be predefined (e.g., the friction factor). These assumptions form a rigidity around the application of numerical methods: they do not provide a straightforward path to adapt to changing system conditions. The first two disadvantages preclude numerical methods as a basis for the control of engineering systems. A majority of model-based control algorithms assume a reduced order or surrogate model for control [7, 8].

In this paper, we introduce Physics-informed State-space Neural Networks for transport phenomena. Physics-Informed Neural Networks (PINNs) are a type of artificial neural network (ANN) used to solve partial differential equations (PDEs) [9]. PINNs combine the expressiveness and flexibility of ANNs with the prior knowledge of physical laws and constraints encoded in PDEs. In a PINN approach, the neural network is trained to approximate the solution of PDEs, given appropriate boundary and initial conditions, by minimizing the difference between the network outputs and the PDE constraints. Thus, PINNs could be a candidate for the solution of a system of PDEs representing transport phenomena. PINNs have several advantages. _Physical priors_: PINNs incorporate prior knowledge about the system through physical constraints. In the context of transport phenomena, prior knowledge is the conservation of mass, momentum, and energy in the spatiotemporal domain. _Mesh-free_: PINNs do not require a mesh or a grid to discretize the spatial domain, as numerical methods do. Instead, PINNs use a neural network to represent the solution in a continuous fashion, without the need for discretization. First, this aspect eliminates the need for mesh generation, which can be time-consuming and challenging. Second, it eliminates discretization error (described in the previous paragraph) arising from the mesh. _Flexible training data_: PINNs can be trained with noisy or scarce data. Thus, sensor measurements, whose outputs have aleatoric sources of noise, may be used to train PINNs. Additionally, the scarcity of data is an essential aspect as the number of sensor locations (i.e., cost) required to utilize PINNs is low.

Specifically, we introduce physics-informed _state-space_ neural networks for transport phenomena. A plurality of previous PINN work focuses on solutions of one or more PDEs, with fixed boundary and initial conditions [10],

\[x\left(z,t\right)=\mathcal{F}\left(z,t\right)\, \tag{1}\]

where \(\mathcal{F}\) is a PINN, \(z\) is a position vector, \(t\) is time, and \(x\) is the field variable(s). However, in a control or state estimation setting, a form representing,

\[x\left(z,t,v\right)=\mathcal{F}\left(z,t,v\right)\, \tag{2}\]

is needed, where \(v\) represents a vector of arbitrary control actions. The control actions modify the boundary conditions (BCs) of the problem. In this setting, the PINN must be able to generalize to a much larger input space. Arnold and King [11] attempted a similar task for a system governed by a single PDE but did not train a monolithic PINN for varying \(v\) - rather multiple PINNs corresponding to a grid of anticipated control actions.

The objective of this work is to develop and demonstrate a method that generates a state-space model by fusing physical measurements from sensors with physical knowledge of transport conservation laws. The paper is organized as follows. In Section 2, the definition and solution of the transport equations, derivation of our proposed method, and its training procedure are presented. Two separate experimental settings were explored. In Section 3, results for modeling a heated channel geometry are presented, and a nonlinear supervisory control method is demonstrated. In Section 4, results for modeling a cooling system are presented, and a demonstration of a capability to use our proposed method for physics-based diagnostics during a system fault. Finally, Section 5 contains a discussion and conclusions on advantages, disadvantages, and additional capabilities to be explored.
## 2 Methods

### Modeling fluid transport

There is a large corpus of equations that model fluid transport. This work focuses on one-dimensional single-phase flow at low-speed, i.e., conditions ubiquitous in engineering systems described in the introduction. Under these conditions, the conservation of mass, momentum, and energy are modeled by the following set of PDEs:

\[\frac{\partial\rho}{\partial t}+\frac{\partial\left(\rho u\right)}{ \partial z} =0\, \tag{3}\] \[\rho\frac{\partial u}{\partial t}+\rho u\frac{\partial u}{ \partial z} =-\frac{\partial p}{\partial z}+\rho g-\frac{f}{D_{h}}\frac{\rho u|u|}{2}\,\] (4) \[\rho C_{p}\frac{\partial T}{\partial t}+\rho C_{p}u\frac{\partial T }{\partial z} =q^{\prime\prime\prime}. \tag{5}\]

In all equations, \(z\) is the spatial dimension, and \(t\) is the temporal dimension. This work uses the field variable set \((p,u,T)\), where \(p\) is the pressure, \(u\) is the one-dimensional velocity, and \(T\) is the temperature. Additionally, there are several parameters. These parameters are either constants or modeled via closure relationships. A closure relationship models density, \(\rho=\rho(p,T)\), pre-defined via experiments. Specific heat capacity is modeled by closure relationships, \(C_{p}=C_{p}(p,T)\), pre-defined via experiments. The friction factor, \(f\), can be obtained via closure relationships (Moody chart) or determined via experiments. The gravitational coefficient, \(g\), is a constant that will depend on the orientation of the component(s) with the gravitational field (and its magnitude). The hydraulic diameter, \(D_{h}\), can be obtained via manufacturer specifications or component measurement. The volumetric heat source, \(q^{\prime\prime\prime}\), depends on the configuration of the thermal energy source, e.g., heater power magnitude and its location.

### Numerical solution

The approach to obtaining numerical solutions to the transport equations is discussed. A solver, System Analysis Module (SAM) [12], is utilized. The SAM code specializes in the high-fidelity simulation of advanced Nuclear Power Plants. As part of its capabilities, it can solve the transport equations, Eqs. (3) to (5). The Finite Element Method (FEM) is utilized to solve the PDEs. FEM works by approximating a continuous solution over a region of interest as a combination of simple functions (e.g., polynomials) that are defined over a finite set of points. The solution is then represented as a set of coefficients describing the simple functions' behavior over the entire region. SAM uses the continuous Galerkin FEM and Lagrange shape functions. This work uses first-order elements using piece-wise linear shape functions with the trapezoidal rule for numerical integration. The Jacobian-Free Newton Krylov method is used to solve the system of equations. The second-order backward difference formula is used for time integration. Additional details on the numerical scheme are provided in [12, SS7].

### PINN solution

PINNs are a type of ANN that can also solve PDEs. The modern PINN concept was initially introduced by Raissi et al. [9] and has since received significant adoption by the scientific machine-learning community - reviewed exhaustively in [10]. PINNs solve parameterized and nonlinear PDEs in the form,

\[\frac{\partial x}{\partial t}+\mathcal{N}(x;\lambda)=0\, \tag{6}\]

where \(x=x(z,t)\) is the solution, and \(\mathcal{N}\) is a nonlinear operator parameterized by \(\lambda\). Each transport equation, Eqs. (3) to (5) can be written in this form. A neural network is formulated to approximate the solution,

\[x(z,t)\approx\mathcal{F}(z,t;\theta),\ z\in\Omega,\ t\in\tau\, \tag{7}\]

where \(\mathcal{F}\) is a neural network that is parameterized by \(\theta\). The solution is box constrained in an n-dimensional space domain, \(\Omega\), and in a one-dimensional time domain, \(\tau\). The space domain specification would depend on the geometry of the transport path and whether multiple dimensions are needed. The specification of the time domain depends on the timescale at which the system reaches a steady state following changes in BCs. For instance, if the system settles slowly, a large \(\tau\) is required, whereas a small \(\tau\) is needed if the system reaches a steady state quickly.

In a purely data-driven approach, a series of measurements would be collected and stored in sets,

\[\mathcal{X}_{m} =\{(z_{i},t_{i})_{i=1}^{N}\}\, \tag{8}\] \[\mathcal{Y} =\{(x_{i})_{i=1}^{N}\}\, \tag{9}\]
where \(z_{i}\) is the position, \(t_{i}\) is the time, and \(x_{i}\) is the field value for the \(i^{\text{th}}\) data point. The set \(\mathcal{X}_{m}\) contains positions and times, where subscript \(m\) denotes measurement data. The set \(\mathcal{Y}\) contains measured field values. The superscript \(N\) is the dataset size. In a supervised learning setting, \(\mathcal{X}_{m}\) is the input dataset, and \(\mathcal{Y}\) is the output dataset. Formally, \(x_{i}\in\mathbb{R}^{m}\), where \(m\) is the number of fields, \(z_{i}\in\mathbb{R}^{n}\cap\Omega\) where \(n\) is the number of spatial dimensions, and \(t_{i}\in\mathbb{R}\cap\tau\) is a scalar. A corresponding loss function is defined,

\[\mathcal{L}_{m}=f_{\text{loss}}\left(\mathcal{F}(\mathcal{X}_{m}),\mathcal{Y }\right), \tag{10}\]

where \(f_{\text{loss}}\) is a loss function (e.g., the mean squared error), and the \(\mathcal{L}_{m}\) denotes the measurement loss. The first argument of the loss function is the predicted value, and the second argument is the target value. Thus, \(\mathcal{L}_{m}\) quantifies the error (or loss) between predictions from the neural network and the measured field values. The loss is then backpropagated to the neural network's parameters, \(\theta\). Each parameter will have a gradient of its value with respect to \(\mathcal{L}_{m}\), known as the local loss gradient. Then, an optimizer will adjust each parameter according to its local loss gradient i.e., change the value of the parameter to decrease the entire network's loss averaged over the current batch. In the purely data-driven setting, the values of \(\theta\) are strictly determined by the measured field values.

In a PINN setting, an additional source of information is added. The dynamics of the system being modeled are approximately known. For example, if a liquid is advected in a network of pipes, the transport equations can inform the neural network that, e.g., energy conservation must hold. Thus, the transport PDEs inform the neural network if there are inconsistencies in its solution. How is this achieved? First, the exact solution in Eq. (6) is replaced by that approximated by the neural network in Eq. (7),

\[\frac{\partial\mathcal{F}}{\partial t}+\mathcal{N}\left(\mathcal{F};\lambda \right)=0\;. \tag{11}\]

To evaluate the PDEs, e.g., Eqs. (3) to (5), the gradients of \(\partial_{t}u,\partial_{z}u\), etc. are needed. In the numerical approach, the gradients are obtained by discretizing the space-time domain and approximating them using finite difference. The PINN approach uses a mesh-free method to evaluate the gradients. This is achieved through _forward mode_ automatic differentiation, where the derivative of any output of \(\mathcal{F}\) is obtainable with respect to any input. The physics-informed loss function is defined as,

\[\mathcal{L}_{p}=f_{\text{loss}}\left(\frac{\partial\mathcal{F}}{\partial t}( \mathcal{X}_{p})+\mathcal{N}\left(\mathcal{F}(\mathcal{X}_{p}),\lambda\right),\emptyset\right)\;, \tag{12}\]

where \(\emptyset\) is a zero set. The target output is always zero to enforce that the neural network satisfies the right-hand side of Eq. (11). It is important to emphasize that the target output is a zero set rather than measured field values. Thus, \(\mathcal{L}_{p}\) can be evaluated at arbitrary inputs sampled from the space-time domain, \((z_{p},t_{p})\sim\Omega\times\tau\). In this way, the known physics of the system regularizes the neural network to the space-time domain outside of the measurement set, \(\mathcal{X}_{m}\). The additional set of inputs is referred to as collocation points,

\[\mathcal{X}_{p}=\left\{(z_{p},t_{p})_{i=1}^{N}\right\}\,. \tag{13}\]

Now, during training, the parameters of the neural network, \(\theta\), can be updated to minimize both \(\mathcal{L}_{m}\) and \(\mathcal{L}_{p}\). Thus, knowledge gathered from measurements of the system is complemented with a priori knowledge of the physics of the system. The resulting neural network is defined as a PINN. Measurements of a facility require the installation of sensors (e.g., thermocouples and pressure transducers) and associated operation and maintenance costs. The amount of sensors and their locations is itself an optimization task. The PINN approach is particularly powerful when system measurements are sparse or physically precluded (e.g., extreme temperature and pressure environments).

### State-space formulation (PSM)

The objective of this work is to apply aspects of PINN to obtain a nonlinear dynamical model as an alternative to traditional linear state-space models. A discrete-time state-space model is represented by,

\[x_{k+1} =f\left(x_{k},v_{k}\right)\;, \tag{14}\] \[y_{k} =g\left(x_{k},v_{k}\right)\;, \tag{15}\]

where \(x_{k}\), \(v_{k}\), and \(y_{k}\) are vectors representing states, inputs, and output measurements at time \(k\). Additionally, \(f\) is the dynamics model, and \(g\) is the measurement model. Linear state-space models are commonly used in control theory because they provide a tractable and simplified approach to analyzing and designing control systems. A linear state-space model utilizes matrices to approximate the \(f\) and \(g\) functions,

\[x_{k+1} =Ax_{k}+Bv_{k}\;, \tag{16}\] \[y_{k} =Cx_{k}+Dv_{k}\;, \tag{17}\]
where \(A\) and \(B\) matrices are the state and input matrices, and \(C\) and \(D\) are the output and control pass-through matrices. The linear state-space representation allows the application of powerful mathematical techniques such as eigenvalue analysis, frequency response analysis, and convex optimization algorithms (e.g., Quadratic Programs for Model Predictive Control). However, many physical systems exhibit nonlinear behavior, including in systems with time-varying parameters or when the system is operating in certain conditions (within steep saddle points). Thus, using linear models can lead to inaccuracies or instability in controlling a nonlinear system.

A nonlinear dynamics model could address the deficiencies of a linear representation. Specifically, a neural network could be used as the nonlinear function representing \(f\) and \(g\) in Eqs. (14) and (15). The previous section presented PINNs as an approach to solving PDEs, Eq. (11). If the PDEs are a good approximation of the dynamics of a physical system, the PINN will approximate the system's state in the form Eq. (1). However, several issues arise. First, the output of the PINN is a continuous representation of the system state \(x=x(z,t)\). A discrete representation could be defined by,

\[x_{k}=\left[\mathcal{F}\left(z_{0},t_{k}\right),\mathcal{F}\left(z_{1},t_{k} \right),...,\mathcal{F}\left(z_{M},t_{k}\right)\right]^{\intercal}\;, \tag{18}\]

where \(z_{0},...,z_{M}\) represent the physical locations at which the system's state is defined, and \(t_{k}\) is the discrete time. Without loss of generality, we can assume that \(t_{k}=k\Delta t\), where \(k\) is an integer indexing time, and \(\Delta t\) is a constant time-step size. Therefore, the propagation of the state would be determined by,

\[x_{0} =\left[\mathcal{F}\left(z_{0},0\right),\mathcal{F}\left(z_{1},0 \right),...,\mathcal{F}\left(z_{M},0\right)\right]^{\intercal}\;, \tag{19}\] \[x_{1} =\left[\mathcal{F}\left(z_{0},t_{1}\right),\mathcal{F}\left(z_{1 },t_{1}\right),...,\mathcal{F}\left(z_{M},t_{k}\right)\right]^{\intercal}\;,\] (20) \[\vdots\] \[x_{N} =\left[\mathcal{F}\left(z_{0},t_{N}\right),\mathcal{F}\left(z_{1 },t_{N}\right),...,\mathcal{F}\left(z_{M},t_{N}\right)\right]^{\intercal}\;. \tag{21}\]

While useful for systems without any inputs, additional modifications are needed to accommodate systems with varying control inputs, i.e., Eq. (2),

\[x_{k}=\left[\mathcal{F}\left(z_{0},t_{k},v_{k}\right),\mathcal{F}\left(z_{1}, t_{k},v_{k}\right),...,\mathcal{F}\left(z_{M},t_{k},v_{k}\right)\right]^{ \intercal}\;, \tag{22}\]

where \(v_{k}\) is a vector representing control inputs at step \(k\). The control inputs can manipulate BCs, e.g., the temperature or velocity at a pipe inlet, or modify jump conditions between locations, e.g., pressure jump due to the presence of a pump. A problem with Eq. (22) is the assumption that the system's initial condition, \(x(z,t=0)\), must always be the same. For a dynamical system, an arbitrary final input value, \(v_{k\rightarrow\infty}\), would result in a new settling point, \(x_{k\rightarrow\infty}\). If we want the state-space model to generalize, a form propagating forward from arbitrary \(x_{0}\) is needed. Therefore, the initial condition itself should be an input for the PINN,

\[x_{k}=\left[\mathcal{F}\left(z_{0},t_{k},x_{0},v_{k}\right),\mathcal{F}\left(z _{1},t_{k},x_{0},v_{k}\right),...,\mathcal{F}\left(z_{M},t_{k},x_{0},v_{k} \right)\right]^{\intercal}\;, \tag{23}\]

where \(x_{0}\) is the initial condition.

When a dynamical system is subject to an initial condition, the system's behavior in the short term is dominated by the initial state. However, as time passes, \(k\gg 1\), the impact of the initial condition on the system's current state becomes progressively smaller. Therefore, a final modification to the formulation is needed,

\[x_{k}=\left[\mathcal{F}\left(z_{0},t_{k},x_{0,k},v_{k}\right),\mathcal{F} \left(z_{1},t_{k},x_{0,k},v_{k}\right),...,\mathcal{F}\left(z_{M},t_{k},x_{0, k},v_{k}\right)\right]^{\intercal}\;, \tag{24}\]

where \(x_{0_{k}}\) is the initial condition of the current time step. Before Eq. (22), it was assumed that \(t_{k}=k\Delta t\). When this assumption is applied to Eq. (24), over a collection of multiple experiments, the value of \(t_{k}\) would arbitrarily correspond to combinations of the inputs, \(\left(x_{0,k},v_{k}\right)\), and outputs, \(x_{k}\). From the perspective of training the neural net \(\mathcal{F}\), there would be no relationship between \(x_{k}\) and \(t_{k}\). However, \(t_{k}\) is still an important variable to estimate the time derivatives needed for physics informing. The timescale at which \(t_{k}\) does have a relationship between the inputs at outputs is the discrete time step, \(\Delta t\). The data would be separated into an initial condition at \(t_{k}=0\), and the next step, \(t_{k}=\Delta t\),

\[x_{0,k} =\left[\mathcal{F}\left(z_{0},0,x_{0,k},v_{k}\right),\mathcal{F} \left(z_{1},0,x_{0,k},v_{k}\right),...,\mathcal{F}\left(z_{M},0,x_{0,k},v_{k} \right)\right]^{\intercal},\;\mathrm{and}, \tag{25}\] \[x_{k} =\left[\mathcal{F}\left(z_{0},\Delta t,x_{0,k},v_{k}\right), \mathcal{F}\left(z_{1},\Delta t,x_{0,k},v_{k}\right),...,\mathcal{F}\left(z_{M},\Delta t,x_{0,k},v_{k}\right)\right]^{\intercal}\;. \tag{26}\]

The relationships between the inputs, \(\left(x_{0,k},v_{k}\right)\), and outputs, \(x_{k}\), in Eq. (24) is visualized in Fig. 1.

The resulting neural network is,

\[x\left(z_{i},t_{i},v_{i},x_{0,i}\right)=\mathcal{F}\left(z_{i},t_{i},v_{i},x_{0,i}\right)\;, \tag{27}\]

predicting each component of Eqs. (25) and (26). To train this neural network, the datasets in Eqs. (8) and (9) are modified,

\[\mathcal{X} =\left\{(z_{i},t_{i},v_{i},x_{0,i})_{i=1}^{N}\right\}, \tag{28}\] \[\mathcal{Y} =\left\{x(z_{i},t_{i},v_{i},x_{0,i})_{i=1}^{N}\right\}, \tag{29}\]
where \(v_{i}\) is the control input, and \(x_{0,i}\) is the initial condition. Formally, \(v_{i}\in\mathbb{R}^{p}\) where \(p\) is the number of control inputs, and \(x_{0,i}\in\mathbb{R}^{q}\) where \(q\) is the number of states modeled. Using this formulation, the input space of the neural network increases by \(\mathcal{O}(10^{p+q})\). However, this step is essential to allow network usage without relying on a heuristic approach, e.g., a tabular database of multiple models where output is obtained by interpolating between \((v_{i},x_{0,i})\). Besides, this would circumvent a desire to show that the network can generalize well. The trained network is referred to as a Physics-informed State-space neural network Model (PSM).

### PSM architecture and training methodology

To train a PSM, an ANN must be constructed. The ANN will represent the operator, \(\mathcal{F}\), in Eqs. (10), (12) and (24). In past research, various PINN architectures have been explored, including Multi-Layer Perceptron (MLP), Convolutional Neural Networks (CNNs), and Recurrent Neural Networks (RNNs) [10]. For this study, we adopt an MLP architecture. The intuition behind this choice is explained by contrasting an MLP's properties against other architectures. Transport system dynamics are inherently Markovian, meaning that the state \(x_{k+1}\) depends solely on the previous state, \(x_{k}\), rather than a sequence of past states. This characteristic makes a state-less architecture, such as an MLP, more suitable than an RNN for this application. Additionally, CNNs are typically employed for tasks involving image or signal processing, where the input exhibits a grid-like structure and exhibits strong spatiotemporal correlations, such as multi-dimensional data. However, in this study, the PSMs are designed to model one-dimensional transport, making an MLP architecture adequate for the task.

Each \(i^{\text{th}}\) layer of the MLP provides an output,

\[y_{j}^{i}=\sigma^{i}\left(w_{jk}^{i}y_{k}^{i-1}+b_{j}^{i}\right)\, \tag{30}\]

where \(w_{jk}^{i}\) is the weight matrix, \(b_{j}^{i}\) is a bias vector, and \(\sigma^{i}\) is the activation function of the \(i^{\text{th}}\) layer. The vectors \(y_{j}^{i}\) and \(y_{k}^{i-1}\) are outputs of the \(i^{\text{th}}\) and \((i-1)^{\text{th}}\) layers. The sizes of the matrices and vectors (denoted by \(j,\ k\)) will depend on the layer-wise size of \(y\). In this work, the ANN is built and trained with PyTorch[13]. The architecture of the neural network and the training workflow is presented in Fig. 2. The current configuration of the neural network is described next. The activation function between all intermediate layers is the hyperbolic tangent. It was chosen because it is nonlinear, is smooth, and has smooth derivatives. After the input layer, there are three "head" layers of size \(s_{H}\), followed by an "intermediate" layer of size \(s_{I}\). The output of the intermediate layer is sent to three separate "tail" layers of size \(s_{T}\). There are three tails for each field variable: pressure, velocity, and temperature. The intuition behind this architecture is that the head layers encode the initial condition and control input into a form used by the tail layers. The tail layers then decode the information separately for each field variable.

As discussed previously, two data sets are used to train the PSM. Measurement data from locations \(z_{i}\in\mathcal{X}_{m}\), Eq. (8), and physics-informing from collocation locations \(z_{i}\in\mathcal{X}_{p}\), Eq. (13). In the context of transport systems, the measurement locations are those where physical measurements of fluid flow exist, e.g., thermocouples or pressure transducers. On the other hand, collocation locations are arbitrarily chosen across the entire physical domain.

#### 2.5.1 Loss formulation

A discussion is presented to explain the differences in the "Measurement Loss", Eq. (10), and "Physics-informed Loss", Eq. (12), in Fig. 2. In conventional PINN formulation, the objective is to solve for the spatiotemporal distribution of field(s) within the boundaries of a domain under the assumptions of a fixed initial condition and BC(s). In general, it is anticipated that there are no measurements inside the spatial domain. However, in transport systems, the direct

Figure 1: A visualization of the nomenclature and their relationships. The red lines indicate variables needed at each step \(k\), i.e., the initial condition and input, \(x_{0,k}\), and \(v_{k}\). The blue arrows indicate the prediction of the next state, \(x_{k}\), by the numerical solver or PSM. The dashed arrows indicate the aliasing of \(x_{k}\) as \(x_{0,k+1}\).
measurements of pressure, temperature, and indirect measurement of velocity, are available through sensors (the type of sensor is dictated by the application, accuracy needed, and fluid properties) at multiple locations along the fluid transport path. Thus, in the PSM formulation, the spatiotemporal distribution of field(s) _between_ the measurement locations must be solved. In this setting, the measurement losses enforce initial conditions in Eq. (25) but also act as BCs in Eq. (26). On the other hand, the physics-informed losses are calculated at random locations. The losses from these locations tune the coefficients of the PSM such that the PDE equations defined in Eq. (11) are satisfied. To tune the coefficients of the PSM, it is important to choose an appropriate loss function.

To train the network, the loss function, \(f_{\mathrm{loss}}\) in Eqs. (10) and (12), must be defined. The loss function quantifies the discrepancy between a target value and a model's prediction. Regression problems in machine learning commonly use the Mean Squared Error (MSE) as the loss function. However, the Log-Cosh loss function, Eq. (31), is utilized in this work. In both equations, \(y_{i}\) is the predicted value, \(\hat{y}_{i}\) is the target, and \(N\) is the total number of samples in a batch.

\[f_{\mathrm{loss}}\left(y_{i},\hat{y}_{i}\right)=\frac{1}{N}\sum_{i=1}^{N}\log \left(\cosh\left(y_{i}-\hat{y}_{i}\right)\right)\;, \tag{31}\]

The Log-Cosh loss function provides benefits over Mean Squared Error (MSE) by being less sensitive to outliers, having bounded gradients for numerical stability, ensuring more reliable convergence in optimization, and performing better with heavy-tailed noise distributions, allowing for quicker, more effective learning and accurate real-world data fitting.

Now that the definitions of the measurement loss, the physics-informed loss, and the loss function have been presented, a discussion on the weighting of the losses follows. The total loss forwarded to the optimizer per epoch is defined as,

\[\mathcal{L}_{\Sigma}=\alpha\mathcal{L}_{m}+\beta\mathcal{L}_{p}\;\big{|}\; \alpha+\beta=1\;, \tag{32}\]

where \(\alpha\) and \(\beta\) are coefficients that weigh the importance of the measurement and physics-informed losses, respectively. The correct settings for \((\alpha,\beta)\) will depend on the degree of certainty in the measurements (assessable from the manufactured accuracy of sensors) and the confidence in the suitability of the PDEs' representation of the system's dynamics. In this work, during training of the PSMs, we set \(\alpha=0.5\) and \(\beta=0.5\), equally weighting both losses. However, further consideration is warranted when applying PSMs to a physical system.

Figure 2: PSM architecture and training workflow. Top: The neural network architecture displaying relative sizes of the MLP layers. Bottom: Measurement and physics-informed losses use output from \(\mathcal{F}\left(\mathcal{X}_{m}\right)\) and \(\mathcal{F}\left(\mathcal{X}_{p}\right)\), respectively. The combined losses are used to update the parameters of the neural network, \(\theta\).
#### 2.5.2 Scaling

Scaling the input and output of a neural network is crucial when processing physical measurements of different units, scales, and ranges, such as pressure, temperature, and velocity. This process facilitates the optimization of the neural network, enabling the model to converge more efficiently and quickly by ensuring that the gradients have a consistent magnitude across all input dimensions. Unscaled inputs with disparate magnitudes can lead to elongated, narrow contours in the loss function landscape, resulting in slow and unstable convergence. In this work, min-max scaling is used,

\[z^{*} =z/z_{\max}\;, \tag{33}\] \[t^{*} =t/t_{\max}\;,\] (34) \[p^{*} =(p-p_{\min})/(p_{\max}-p_{\min})\;,\] (35) \[u^{*} =(u-u_{\min})/(u_{\max}-u_{\min})\;,\] (36) \[T^{*} =(T-T_{\min})/(T_{\max}-T_{\min})\;,\] (37) \[\rho^{*} =(\rho-\rho_{\min})/(\rho_{\max}-\rho_{\min})\;, \tag{38}\]

for transforming position, time, pressure, velocity, temperature, and density, respectively. The value for \(z_{\max}\) equals the maximum path length for the liquid. The value for \(t_{\max}\) equals the discrete time-step between measurements, \(\Delta t\) in Eq. (26). The minimum and maximum values for the pressure, velocity, temperature, and density are based on the ranges of data that are collected from experiments.

#### 2.5.3 Training procedure

The PSM training procedure is listed in Alg. 1. During each epoch, the weights are updated from a fixed measurement dataset. Whereas the collocation points are updated at each epoch. A step learning rate scheduler was used, with a rate decay of 0.5 every 50 epochs. All models were trained for 500 epochs with a batch size of 2048.

```
1:Input: Measurement data \((\mathcal{X}_{m},\mathcal{Y})\), Governing equations \(G\), Loss coefficients \((\alpha,\beta)\), Number of epochs \(N\), Batch size \(B\), Learning rate scheduler \(\mathcal{S}(\cdot)\)
2:Initialize: Network \(\mathcal{F}\) with random weights \(\theta\)
3:for\(n=1,2,\ldots,N\)do
4: Update learning rate \(\eta\leftarrow\mathcal{S}(n)\)
5:for each mini-batch \((\mathcal{X}_{m,i},\mathcal{Y}_{i})\subset(\mathcal{X}_{m},\mathcal{Y})\) with size \(B\)do
6: Compute predictions \(\hat{\mathcal{Y}}_{i}=\mathcal{F}(\mathcal{X}_{m,i})\)
7: Compute measurement loss \(\mathcal{L}_{m}(\hat{\mathcal{Y}}_{i},\mathcal{Y}_{i})\) using Eq. (10)
8: Randomly sample \(B\) collocation points \(\mathcal{X}_{p,i}\), and extract \((x_{0,i},v_{i})\) from \(\mathcal{X}_{m,i}\)
9: Compute the residual values \(\mathcal{R}_{i}=G(\mathcal{F},\mathcal{X}_{p,i},x_{0,i},v_{i})\)
10: Compute physics-informed loss \(\mathcal{L}_{p}(\mathcal{R}_{i},\emptyset)\) using Eq. (12)
11: Compute total loss \(\mathcal{L}_{\Sigma}=\alpha\mathcal{L}_{m}+\beta\mathcal{L}_{p}\)
12: Update network weights using optimizer: \(\theta\leftarrow\theta-\eta\nabla_{\theta}L_{\Sigma}(\theta)\)
13:endfor
14:endfor
15:return Trained PSM model \(\mathcal{F}\)
```

**Algorithm 1** Training a PSM

## 3 Results: applying PSMs for transport in a Heated Channel

Two experiments are designed to investigate PSMs' performance for transport phenomena. The first experiment is a system of three conduits in series (referred to as pipes). The first and last pipes will be adiabatic and \(1.0\,\mathrm{m}\) long. The middle pipe will be heated with a constant volumetric heat source of \(q^{\prime\prime\prime}=50.0\,\mathrm{MW}/\mathrm{m}^{3}\), and \(0.8\,\mathrm{m}\) long. All pipes have a total cross sectional area of \(0.449\,\mathrm{m}^{2}\) and a hydraulic diameter of \(2.972\,\mathrm{mm}\). No solids are modeled - only the internal flow area. Each pipe is discretized into ten elements. The inlet of the first pipe will have varying velocity and temperature BCs,

\[u(z =z_{\min},t) =u_{\mathrm{in}}(t)\;, \tag{39}\] \[T(z =z_{\min},t) =T_{\mathrm{in}}(t)\;. \tag{40}\]
The trajectory of the BCs is defined by linear ramps that vary in ramp rate, maximum/minimum amplitude, and rest time between manipulations. These parameters are varied to encourage regularization of the PSM. Additionally, the outlet of the third pipe will be a constant pressure BC,

\[p(z=z_{\mathrm{max}},t)=p_{\mathrm{out}}. \tag{41}\]

In this work, the pressure BCs are presented as absolute pressure; however, the pressure field is presented as gage pressure. The experiment configuration is visualized in Fig. 3. A visualization of the evolution of the field variables is presented in Fig. 4. The PSM will be tasked with predicting these spatiotemporal distributions by learning from a few measurement locations and the approximate PDEs of the system. Hereafter, this experiment will be referred to as the "heated channel" configuration.

To close the transport equations, Eqs. (3) to (5), the thermophysical properties of the fluid must be defined. This work utilizes a molten salt, \(\mathrm{LiF}-\mathrm{BeF}_{2}\), known colloquially as 'flibe.' The reason is that flibe is the leading candidate for the primary working fluid of next-generation molten salt reactor designs currently under development [14]. The density and specific heat capacity of flibe are defined by [15],

\[\rho=2413-0.488T\;\mathrm{[kg/m^{3}]}\;, \tag{42}\] \[C_{p}=2414\;\mathrm{[J/(kg\,K)]}\;, \tag{43}\]

where \(T\) is the fluid temperature in Kelvin. These correlations are adopted to calculate the physics-informed loss in Eq. (12).

Figure 4: An overview of the training dataset for the heated channel configuration. Left: Velocity and temperature BCs at the inlet that are manipulated in time. The experiment chosen to display the solution is annotated (”Exp. 8”). Right: Numerical solution of all field variables over the entire spatiotemporal domain.

Figure 3: Configuration of the pipes in series. Left: Layout of the pipes and identifications of the locations of imposed BCs and heat source term. Right: A list of the fixed outlet pressure BC, and ranges for the inlet BCs.
### Evaluation with test data

The results for the heated channel are discussed in this section. A sample of the control input permutations and numerical solution of the field variables were presented in Fig. 4. In this setting, the PSM must accurately propagate the system's state, \(x_{k}\), from a given initial condition, \(x_{0,k=0}\), and input trajectory, \(v_{k=0,..N}\). This problem is a classical advection problem where the inlet boundary conditions are transported through the pipe incur: pressure losses due to friction from the pipe, modifications in mass flow rate (velocity) due to boundary conditions, and changes in temperature due to energy deposition.

In Fig. 5, a rollout of a single input, selected from the _test_ dataset, i.e., data that is not used to train the models, is presented. Subplots A and B show the input trajectory of the velocity and temperature at the inlet, \(u_{\mathrm{in}}\) and \(T_{\mathrm{in}}\). Subplots C-K show pressure, velocity, and temperature at fixed positions vs. time. Subplots L-T show the fields at fixed times vs. position. In subplots C-T, results from the numerical solver ("Num. Sol."), the PSM model, and an Artificial Neural Network (ANN) model are displayed. The PSM model is trained with Eq. (32), which contains both measurement and physics-informed loss, Eqs. (10) and (12), respectively. Whereas the ANN model is trained only with the measurement loss. Thus, the ANN results represent the classical, purely data-driven approach to modeling. In contrast, the PSM results combine knowledge from data with a priori knowledge of the system's dynamics in the form of transport PDEs.

The ANN appears to perform on par with the PSM. Comparing the pressure and velocity spatiotemporal predictions, indeed, the ANN does perform as well as the PSM. However, comparing the temperature distributions, there is a significant error introduced by the ANN. The errors occur at \(1.0\leq z\leq 1.8\,\mathrm{m}\), i.e., the heated pipe's location. This observation is more explicitly visualized in Fig. 6, which compares spatiotemporal error in predicting temperature averaged over all test datasets. The PSM significantly outperforms the ANN in the prediction of the temperature field.

Why does the ANN fail to accurately predict the temperature field only? This question can be answered by understanding the underlying physics and the sensor measurement locations.

The advection of pressure and velocity depends on the frictional losses and the inlet mass flow rate. Because the friction factor and pipe geometry are constant, we expect a constant pressure gradient with respect to position. This is displayed in Subplots L-N in Fig. 5. The measurement of the pressure and velocity fields occurs in the first and last unheated pipes. Therefore, the ANN has to simply interpolate between the measurement locations of pressure before and after the heated section to correctly estimate the pressure field, and be able to modify the gradient according to the change in inlet velocity mass flow rate. The latter is provided as an input (\(v_{k}\)) to the ANN. The results show that the ANN model can accomplish these tasks well for the pressure field.

For the velocity field, we expect the spatial distribution to be nearly constant, except after the heated region. There is a change in the density because of the heating (Eq. (42)), which causes the velocity to change in order to conserve mass (Eq. (3)). This causes a small error in ANN prediction, noticeable upon close inspection of Subplots O-Q in Fig. 5. Therefore, the performance for velocity is still acceptable.

Figure 5: Contrasting the performance of PSM and ANN models in predicting the evolution of pressure, velocity, and temperature for the heated channel. Subplots A-B: control input setting from the test dataset. Subplots C-K: field temporal evolution at fixed positions. Subplots L-T: field spatial distribution at fixed times.
However, for the temperature field, there is a large discrepancy between ANN and PSM performance within the heated region. The temperature field is affected by energy deposition (\(q^{\prime\prime\prime}\) in Eq. (5)), and inlet BCs. When \(q^{\prime\prime\prime}=0\), there is no gradient in temperature. When \(q^{\prime\prime\prime}=c\), a constant, there is a linear gradient in temperature. The PSM correctly predicted both situations, as shown in Subplots I-K and R-T in Fig. 5. However, the ANN has no knowledge of \(q^{\prime\prime\prime}\), and thus it interpolates between the sensor readings. Therefore, the ANN model can only correctly predict temperature distributions in regions where \(q^{\prime\prime\prime}=0\), i.e., there is no temperature gradient and any modification is provided directly by the inlet BC.

### Training with noisy measurements

The data from sensor measurements in a physical plant are noisy. The noise stems from inherent randomness during fluid transport & signal measurement itself. It is important to quantify the performance of PSMs in such environments. The introduction of homoscedastic and heteroscedastic noise is studied. The former is introduced by,

\[x_{\mathrm{noisy}}=x+\sigma\cdot\epsilon\;, \tag{44}\]

where \(x_{\mathrm{noisy}}\) is the noisy signal, \(x\) is the noise-free signal, \(\sigma\) is the noise standard deviation, and \(\epsilon\) is a random noise vector sampled from a standard normal distribution. This type of noise models scenarios where the variance is constant across all magnitudes. Heteroscedastic noise is introduced by,

\[x_{\mathrm{noisy}}=x+\sqrt{|x|\cdot k}\cdot\epsilon\;, \tag{45}\]

where \(k\) is the variance factor. The variance factor is a parameter that controls the relationship between the signal's magnitude and the added noise's variance. Thus, this type of noise models scenarios where the noise variance is not constant but depends on the underlying signal's magnitude. The noisy signals are visualized in Fig. 7. During training, the signals are manipulated randomly for each batch. Noise is added to the sensor measurements and the initial conditions (\(x_{k}\) and \(x_{0,k}\), respectively, in Eq. (26)). It is assumed that the system's actuation is not noisy (i.e., the control signals, \(v_{k}\) are noise-free).

The results from introducing noise are discussed next. The training procedure remains the same, except that artificial noise is added during each batch. In Table 1, the results are tabulated for both types of noise, and varying magnitudes of noise factors (\(\sigma\) in Eq. (44), \(k\) in Eq. (45)). First, inspecting the mean RMSE values, the PSM has a significant advantage in predicting velocity and temperature fields. There is a lower advantage for the pressure fields. This might be because the nominal RMSE values for pressure are already very low when contrasted with the pressure range (\(\Delta P\)) across all experiments. Inspecting the maximum RMSE values, the PSM retains a significant advantage in predicting velocity and temperature fields, with some anomalies noted for \(\sigma=0.05\) and \(k=0.05\). The PSM's advantage for the pressure field is relatively moderate, likely because it has a simple spatiotemporal relationship in the current experimental setting. To summarize, the results show that even when noise is added to the signals, PSMs have a significant advantage over a purely data-driven model.

Figure 6: Comparison of root-mean-square error (RMSE) in the prediction of temperature for PSM (left) and ANN (right) models. The RMSE values are averaged over all test datasets. The heated region is present between the dashed red lines.
### PSM Capability: Non-linear Supervisory Control

In the previous section, the PSM architecture was shown to approximate the spatiotemporal distribution of all the fields accurately. Because the PSM is an end-to-end differentiable model, we can calculate the gradients for inputs and outputs using auto differentiation. In this section, the capability of a PSM to provide non-linear supervisory control is demonstrated. The Reference Governor (RG) supervisory control scheme is adopted to demonstrate this capability. The RG algorithm enforces pointwise-in-time constraints on a system by modifying the reference input setpoints [16].

The RG algorithm is summarized next. At each time-step \(k\in\mathbb{Z}_{+}\), the RG receives a reference input, \(r_{k}\in\mathbb{R}^{m}\), where \(m\) is the number of inputs. Depending on the system's current state, expected evolution, and imposed constraints, an admissible input \(v_{k}\in\mathbb{R}^{m}\),

\[v_{k}=v_{k-1}+\kappa_{k}\left(r_{k}-v_{k-1}\right)\;, \tag{46}\]

is obtained. The admissible input is then sent to lower-level controllers of the system (e.g., a PID controller). In the Scalar RG (SRG) formulation, \(\kappa_{k}\in[0,1]\) is a _scalar_ that governs admissible changes to the inputs, such that a complete rejection, \(v_{k}=v_{k-1}\), acceptance, \(v_{k}=r_{k}\), or an intermediate change, \(v_{k-1}\leq v_{k}\leq r_{k}\), is possible.

Next, how the SRG algorithm determines \(\kappa_{k}\) is summarized. The RG uses a linear state-space representation, Eqs. (16) and (17), and then imposes constraints on the output variables, \(y_{k}\in Y\), where a set of linear inequalities defines \(Y\). Then, a construct referred to as the Maximal Output Admissible Set (\(O_{\infty}\)) is defined - it is the set of all \(x_{k}\), and constant

\begin{table}
\begin{tabular}{|l c c c c c c|}  & \(\varepsilon_{P,\mathrm{ANN}}[\mathrm{Pa}]\) & \(\varepsilon_{u,\mathrm{ANN}}[\mathrm{mm/s}]\) & \(\varepsilon_{T,\mathrm{ANN}}[^{\circ}\mathrm{C}]\) & \(\frac{\varepsilon_{P,\mathrm{PSM}}}{\varepsilon_{P,\mathrm{ANN}}}[\%]\) & \(\frac{\varepsilon_{u,\mathrm{PSM}}}{\varepsilon_{u,\mathrm{ANN}}}[\%]\) & \(\frac{\varepsilon_{T,\mathrm{PSM}}}{\varepsilon_{T,\mathrm{ANN}}}[\%]\) \\ \hline \multicolumn{7}{|c|}{Mean} \\ \hline Noise-free & 0.51 & 0.17 & 1.01 & 58 & 52 & 6 \\ Homos. \(\sigma=0.01\) & 1.25 & 0.25 & 1.15 & 12 & 25 & 3 \\ Homos. \(\sigma=0.05\) & 0.46 & 0.22 & 0.98 & 143 & 35 & 17 \\ Homos. \(\sigma=0.10\) & 0.83 & 0.23 & 1.01 & 91 & 51 & 28 \\ Heteros. \(k=0.01\) & 0.71 & 0.25 & 1.15 & 94 & 41 & 24 \\ Heteros. \(k=0.05\) & 0.94 & 0.27 & 1.10 & 94 & 68 & 27 \\ Heteros. \(k=0.10\) & 1.07 & 0.34 & 0.97 & 89 & 63 & 33 \\ \hline \multicolumn{7}{|c|}{Maximum} \\ \hline Noise-free & 4.82 & 1.25 & 5.81 & 113 & 52 & 31 \\ Homos. \(\sigma=0.01\) & 8.97 & 1.24 & 6.08 & 10 & 43 & 13 \\ Homos. \(\sigma=0.05\) & 6.17 & 1.24 & 7.00 & 53 & 104 & 14 \\ Homos. \(\sigma=0.10\) & 6.49 & 2.38 & 6.87 & 78 & 35 & 27 \\ Heteros. \(k=0.01\) & 4.80 & 1.34 & 7.70 & 81 & 51 & 19 \\ Heteros. \(k=0.05\) & 5.47 & 1.42 & 6.54 & 95 & 120 & 35 \\ Heteros. \(k=0.10\) & 6.31 & 1.97 & 5.46 & 76 & 69 & 38 \\ \hline \multicolumn{7}{|c|}{Range across all experiments: \(\Delta P=531\,\mathrm{Pa}\), \(\Delta u=200.0\,\mathrm{mm/s}\), \(\Delta T=92.6\,^{\circ}\mathrm{C}\)} \\ \end{tabular}
\end{table}
Table 1: Comparison of the mean and maximum RMSE for field variables using the ANN across all test datasets. The relative RMSE error for PSMs is displayed in the last three columns: the values are color-coded to highlight significant performance improvement when below 75 %, and a degradation in performance if above 100 %. The range of field values over all experiments is displayed at the bottom.

Figure 7: Example of the manipulation of sensor measurements for all fields. During training, the signals are manipulated randomly for each batch.
input \(\tilde{v}\), such that,

\[O_{\infty}=(x_{k},\tilde{v}):y_{t+k}\in Y,v_{t+k}=\tilde{v},\forall k\in\mathbb{Z }^{+}. \tag{47}\]

The \(k\to\infty\) assumption of \(O_{\infty}\) is generally relaxed to a suitably large finite horizon, \(T\). Thus, at each time-step \(k\), the SRG algorithm determines the admissible \(\tilde{v}\), and therefore \(\kappa_{k}\), such that the system remains in \(O_{\infty}\). The \(O_{\infty}\) is constructed by a roll-out of Eqs. (16) and (17), detailed further in [16, 82]. For a fixed set of state-space matrices and constraints, the matrices needed to construct \(O_{\infty}\) can be precomputed and stored for each RG evaluation.

A drawback of the SRG formulation is that \(\kappa_{k}\) is a scalar. If \(v_{k}-r_{k}\) is multi-dimensional and \(\kappa_{k}<1\), input movement in all dimensions is bounded. In this work, we utilize the Command Governor (CG), a variant of the RG which selects \(v_{k}\) by solving the quadratic program,

\[v_{k}= \operatorname*{argmin}_{v_{k}}\left\|v_{k}-r_{k}\right\|_{Q}^{2}\,\] \[\mathrm{s.t.} (x_{k},v_{k}=\tilde{v})\in O_{\infty}\, \tag{48}\]

where \(Q\) is a positive definite matrix signifying the relative importance of each \(v_{k}\) component. The quadratic program is solved by the CVXPY library [17]. Using the CG, constraints can be enforced by manipulating \(v_{k}\) arbitrarily in \(\mathbb{R}^{m}\), constrained by \(O_{\infty}\).

To utilize the CG, a linear representation of the system is needed. This work assumes full-state feedback, \(C=\mathbb{I}\), and no control pass-through, \(D=\mathbf{0}\). Because the PSM is end-to-end differentiable, we can approximate the nonlinear model, Eq. (14), as a linear representation using Jacobian linearization,

\[x_{k+1} \approx x_{00}+A\delta x_{k}+B\delta v_{k}\, \tag{49}\] \[\delta x_{k} =x_{k}-x_{00}\,\] (50) \[\delta v_{k} =v_{k}-v_{00}\, \tag{51}\]

where \(x_{00}\) and \(v_{00}\) are state and input values at the linearization point. The matrices \(A\), and \(B\) are approximated by the Jacobians of \(\mathcal{F}\), Eq. (27), at the linearization point,

\[A =\left[\frac{\partial\mathcal{F}_{x_{k+1}}}{\partial x_{k,0}}, \frac{\partial\mathcal{F}_{x_{k+1}}}{\partial x_{k,1}},...,\frac{\partial \mathcal{F}_{x_{k+1}}}{\partial x_{k,N}}\right]^{\intercal}_{(x_{00},v_{00})}\, \tag{52}\] \[B =\left[\frac{\partial\mathcal{F}_{x_{k+1}}}{\partial v_{k,0}}, \frac{\partial\mathcal{F}_{x_{k+1}}}{\partial v_{k,1}},...,\frac{\partial \mathcal{F}_{x_{k+1}}}{\partial v_{k,M}}\right]^{\intercal}_{(x_{00},v_{00})}\, \tag{53}\]

where \(N\) and \(M\) are the number of states and inputs. In summary, the RG was introduced as the framework to enforce constraints, followed by a specific version, the CG, which can manipulate multiple inputs to enforce constraints on a system. Then, the approach to utilize the PSM to approximate the linear state-space matrices required by the CG was introduced. The algorithm to utilize this procedure for supervisory control is listed in Alg. 2.

```
1:Input: PSM model \(\mathcal{F}\), Requested set point trajectory \(r\), Time-dependent constraints on states \(Y_{k}\), Input matrix \(Q\), Matrices \(A\) and \(B\) update interval \(\gamma\)
2:Initialize: Initial condition of system \(x_{0}\), Iteration counter \(i\gets 0\)
3:for\(k=1\) to \(\text{length}(r)\)do
4:if\(k\mod\gamma=0\)then
5: Obtain matrices \(A\) and \(B\) by Jacobian Linearization of \(\mathcal{F}\) at state \(x_{k},v_{k}\)
6: Update the \(O_{\infty}\) matrices with the current linear state-space model
7:endif
8:if Constraint update then
9: Update the \(O_{\infty}\) matrices with the current \(Y_{k}\)
10:endif
11: Formulate the quadratic program in Eq. (48) with inputs \(r_{k},x_{k}\)
12: Solve the quadratic program to obtain the admissible set point \(v_{k}\)
13:Return Admissible set point \(v_{k}\) to controlled system
14: Update state of system \(x_{k}\gets x_{k+1}\)
15:endfor
```

**Algorithm 2** Sequential Neural Command Governor (NCG) using PSMs

Next, results demonstrating Alg. 2 are presented. In Fig. 8, an input permutation from the test dataset is chosen randomly to manipulate the inlet BCs. The variation of the temperature BC, \(T_{\rm in}\), causes a significant temperature
rise after the heated channel, observable by comparing subplots I-K. It is desirable to constrain excessive increases in temperature to avoid damage to equipment or mitigate corrosion. Thus, a temperature constraint is assigned at a location after the heated channel, \(z=$2.3\,\mathrm{m}$\), shown in subplot K. The NCG algorithm admits input changes to enforce these constraints, shown in subplots A and B. The results demonstrate that the approach is successful in meeting the enforced constraints. There are a few key takeaways. First, the NCG algorithm independently manipulates each input value, minimizing the distance from the requested input trajectories, i.e., Eq. (48). Second, a time-dependent constraint is demonstrated (annotated in subplot K). Accommodating time-dependent constraints enables flexibility during operation - it is envisioned that, as part of a comprehensive autonomous operation framework, an external algorithm could update the constraints enforced on a system depending on its health condition. Lastly, it is important to note that in this demonstration, the PSM is only used to update the \(A\) and \(B\) matrices - the numerical solver is used as the environment with which the NCG algorithm interacts.

## 4 Results: applying PSMs for transport in a Cooling System

In addition to the heated channel configuration, a "cooling system" configuration was adopted. A series of six pipes are arranged in a loop. All pipes have a total cross sectional area of \(0.449\,\mathrm{m}^{2}\) and a hydraulic diameter of \(2.972\,\mathrm{mm}\). No solids are modeled - only the internal flow area. All pipes are \(1.0\,\mathrm{m}\) in length, except for the pipes immediately after the heated and cooled sections, which are \(2.0\,\mathrm{m}\) long. The pipes are discretized into 10 elements per meter. A pump provides the head necessary to drive the fluid flow. In addition to a volumetric heat source, a volumetric heat sink is also present. This configuration is commonly found in advanced cooling systems, such as those used in high-performance electronic devices, nuclear reactors, and aerospace applications. The loop configuration is visualized in Fig. 9. There are two independent control inputs,

\[q_{\mathrm{in}}^{\prime\prime\prime} =q_{\mathrm{in}}^{\prime\prime\prime}(t)\;, \tag{54}\] \[q_{\mathrm{out}}^{\prime\prime\prime} =-q_{\mathrm{in}}^{\prime\prime\prime}(t)\;,\] (55) \[\Delta p_{\mathrm{pump}} =\Delta p_{\mathrm{pump}}(t), \tag{56}\]

where \(q_{\mathrm{in}}^{\prime\prime\prime}\) is the volumetric heat generation rate, \(q_{\mathrm{out}}^{\prime\prime\prime}\) is the volumetric cooling rate, and \(\Delta p_{\mathrm{pump}}\) is the pressure head provided by the pump. A control strategy is adopted that synchronizes the heating and cooling rates to increase/decrease the system's heat generation/cooling capability while maintaining constant temperatures throughout the loop. Both independent control inputs (\(q_{\mathrm{in}}^{\prime\prime\prime}(t),\Delta p_{\mathrm{pump}}(t)\)) are a function of time and will be defined by linear ramps that vary in ramp rate, maximum/minimum amplitude, and rest time between manipulations. Additionally, there is one boundary condition,

\[p(z=z_{\mathrm{set}},t)=p_{\mathrm{loop}}\;, \tag{57}\]

which is a static pressure boundary that defines the relative value of the pressure in the loop.

Figure 8: Demonstration of the NCG algorithm. The PSM instantiates the models needed to utilize the supervisory control algorithm. The numerical solver is used as the environment with which the NCG algorithm interacts. Subplots A-B: Control input permutations. The input value without constraints is the reference input trajectories for the NCG algorithm. Subplots C-K: Field values at fixed locations. Subplot K: A constraint is assigned at this location and manipulated as a function of time.
Unlike the heated channel configuration, the loop configuration is more challenging from a dynamics modeling perspective. First, there is no time-dependent manipulation of the field variables via boundary conditions. The manipulation of the pump head and heat source indirectly impact the field variables. Their evolution will have a spatiotemporal effect dictated by the transport equations, Eqs. (3) to (5). A visualization of the evolution of field variables is presented in Fig. 10. The pump head trajectory strongly impacts the velocity of the fluid. Due to the control strategy, the temperature of the fluid remains approximately constant with respect to time but has a spatial variation.

### Evaluation with test data

The performance of the PSM and ANN models for the cooling system is discussed next. A sample of the control input permutations and the numerical solution of the field variables are presented in Fig. 11. Through visual inspection, it is apparent that the ANN struggles with accurately predicting the qualitative evolution of all the field variables. Because the ANN relies on a purely data-driven approach, the neural network correctly fits the field values about the spatial locations corresponding with sensors. Yet, it fails to correctly interpolate the field values between the sensor locations

Figure 10: An overview of the training dataset for the loop configuration. Left: Pump pressure jump and heat source inputs manipulated in time. The experiment chosen to display the solution is annotated (”Exp. 8”). Right: Numerical solution of all field variables over the entire spatiotemporal domain.

Figure 9: Configuration of pipes in a loop. Left: The imposed pump BC location and heat source/sink terms. Right: A list of the ranges for the pump pressure jump, heat source/sink terms, and loop pressure setting.
because the problem is more complicated than the advection-dominated heated channel. Whereas, the PSM has good quantitative and qualitative spatiotemporal prediction of all field variables.

The contrast between the ANN and PSM models is illustrated further in Fig. 12. The figure presents error in predicting all field variables averaged over all test datasets. The PSM model offers a significant performance advantage. As the prediction horizon for the ANN increases, a significant error is introduced in all fields. However, even for small time horizons, \(t\leq 100\,\mathrm{s}\), there are significant errors in the prediction of temperature and pressure. In summary, the results in this section show that in a complex experimental setting, a purely data-driven approach, i.e., the ANN model, overfits to sensor measurements and results in unphysical predictions over the entire spatiotemporal domain. In contrast, the PSM can generalize well to test data and achieve low error.

Figure 11: Contrasting the performance of PSM and ANN models in predicting the cooling system’s evolution of pressure, velocity, and temperature. Subplots A-B: control input setting from the test dataset. Subplots C-K: field temporal evolution at fixed positions. Subplots L-T: field spatial distribution at fixed times.

Figure 12: Comparison of the RMSE in predicting the spatiotemporal evolution of all field variables, averaged over all test datasets. Subplots A-C display error for the PSM model, and subplots D-F display error for the ANN model.
### PSM Capability: Diagnostics

Diagnostics, in the context of this work, refers to the systematic identification of system degradation or faults. A key benefit of utilizing PSMs is their ability to incorporate the governing physics of the system through PDEs. These PDEs produce residuals - differences between observed and predicted behavior - that can be inspected for diagnostic purposes.

In a fully trained PSM, residuals from the mass, momentum, and energy PDEs, Eqs. (3) to (5), stabilize during the system's nominal fault-free operation. If a component degrades during operation, these residuals will change. An algorithm can then use these changes in the residual signature to identify and classify the nature of the degradation.

To validate this approach, a test case was designed to simulate a degradation in the cooling system's performance. In this case, an abrupt partial blockage of the system occurs in Pipe 3 (immediately before the heat sink). The blockage was simulated by increasing the friction factor, \(f\) in Eq. (4), of that pipe by an order of magnitude. In Alg. 3, a procedure is proposed to detect degradation and use a classification algorithm to diagnose the incident. The algorithm requires continuous collection and analysis of new datasets, \((\mathcal{X}_{\mathrm{new}},\mathcal{Y}_{\mathrm{new}})\), and analyzing them. It operates using two PSMs: a 'nominal model', which represents the fault-free system, and a 'twin model' whose parameters are updated once degradation is detected. The update is triggered by a parameterized, \(\tau\), in Line 6.

After the twin model is updated, residuals from both PSMs are analyzed using a classification model to infer the type of degradation that has occurred. The development of this classification model is outside the scope of the current work.

```
1:Input: Original dataset \((\mathcal{X},\mathcal{Y})\), New dataset \((\mathcal{X}_{\mathrm{new}},\mathcal{Y}_{\mathrm{new}})\), Nominal PSM model \(\mathcal{F}_{\mathrm{nom}}\), Threshold \(\tau\), Classification model \(\mathcal{C}\)
2:Initialize: Twin PSM model \(\mathcal{F}_{\mathrm{deg}}\), Degradation flag latch\(\leftarrow\) False
3:for all sample \((x_{k},y_{k})\) in \((\mathcal{X}_{\mathrm{new}},\mathcal{Y}_{\mathrm{new}})\)do
4: Predict \(\hat{y}_{k}=\mathcal{F}_{\mathrm{nom}}(x_{k})\)
5: Compute \(E=\frac{1}{n}\sum_{k=1}^{n}(\hat{y}_{k}-y_{k})^{2}\)
6:if\(E>\tau\)and latch is Falsethen
7: Set latch to True
8:endif
9:if latch is True then
10: Train \(\mathcal{F}_{\mathrm{deg}}\) on \((x_{k},y_{k})\) without physics-informing
11:endif
12:endfor
13: Compute nominal residuals as \(r=\text{PDEs}(\mathcal{F}_{\mathrm{nom}}(x))\) for \(x\) in \(\mathcal{X}\)
14: Compute degraded residuals as \(r_{\mathrm{deg}}=\text{PDEs}(\mathcal{F}_{\mathrm{deg}}(x))\) for \(x\) in \(\mathcal{X}\)
15: Create combined residual dataset \(\mathcal{R}=\{r,r_{\mathrm{deg}}\}\)
16: Predict degradation type \(d\) using \(\mathcal{C}\) on \(\mathcal{R}\): \(d=\mathcal{C}(\mathcal{R})\)return\(d\)
```

**Algorithm 3** Degradation Detection and Classification

The application of Alg. 3 to the data generated from the degradation test case is discussed next. In Fig. 13, the scaled residual values are contrasted for nominal and degraded models. The residuals are presented as a function of position. Within the cooling system, Pipe 3 is located at \(4.0\,\mathrm{m}\leq z\leq 5.0\,\mathrm{m}\). There are a few important observations to consider. First, the nominal model has a relatively constant spatial distribution for the momentum residual, \(r_{\mathrm{mom}_{\cdot}}\). Whereas, for the degraded model, the spatial distribution shows a qualitative shift in the residual around the location of Pipe 3. Second, the mass and energy residuals do not show a significant shift in the spatial distribution of residuals. From these observations, a classification algorithm can be trained to identify the degradation that has occurred. In summary, in this section, an algorithm was proposed and utilized on a test case to show that PSMs can provide meaningful data to a diagnostics algorithm to identify system degradation or fault.
## 5 Discussion and Concluding Remarks

Autonomous systems are envisioned to achieve real-time optimization, be flexible in operation, and be fault tolerant. These capabilities necessitate a model-based approach. However, purely data-driven methods lack physical constraints like mass conservation. Consequently, model-free approaches, such as PID control, dominate the operation of modern transport-dominated systems such as chemical, biomedical, and power plants. To address this challenge, we propose the adoption of PSMs.

A PSM is a deep neural network trained by fusing sensor data with physics-informing using the components' PDEs. The result is a physics-constrained end-to-end differentiable forward dynamics model. Two separate experiments were designed, in silico, to demonstrate PSM capabilities for transport-dominated systems. The first experiment was a heated channel, and the second was a cooling system loop. The datasets used to train the PSM models are available online1. Evaluation on test data indicates that the PSM approach is quantitatively and qualitatively more accurate than a purely data-driven approach, such as an ANN. Additionally, both homo- and heteroscedastic noise were introduced during training to demonstrate the applicability of PSMs to a physical plant.

Footnote 1: Dataset will be released on journal acceptance.

Because PSMs are an end-to-end differentiable state-space model, there are several use cases:

1. The first demonstration used a PSM model to construct a supervisory controller (Section 3.3). To achieve this, the proposed NCG algorithm, Alg. 2, used a linear state-space representation sequentially updated through Jacobian Linearization of the PSM. The capability was utilized to demonstrate both constant and time-dependent constraint enforcement. Thus, PSMs can be utilized for control schemes that require a non-linear or linear state-space representation.
2. The second demonstration proposed a diagnostics algorithm that utilized the PSM as a basis (Section 4.2). Conservation laws, in the form of PDEs, can be evaluated using the PSM. The proposed algorithm, Alg. 3, uses residuals from each of the PDEs, which are available as a function of location. The demonstration simulated a partial blockage of a pipe in the cooling system. Two separate PSMs were used to obtain PDE residuals, one trained on the nominal dataset and one that had been transfer-learned on the degraded dataset. It was demonstrated that the differences in the residuals of the momentum equation can then be used as inputs to a classification algorithm.
3. An additional use-case of PSMs, not covered in this work, is as a basis for Digital Twins (DTs). DTs require that the "digital" model used to represent the physical asset is updated frequently through feedback from measurements to match drifts from anticipated dynamics e.g., due to wear and tear [18]. Because the PSM is end-to-end differentiable, this feature can be accommodated through online learning of PDE coefficients.

There are potential drawbacks of the PSM that warrant further research. First, this work did not rigorously study the number and placement of sensors (measurement locations). The optimal choice would be highly dependent on the characteristics of the system. If a PSM is applied to an existing facility, adding additional sensors may not be possible. Second, the current work intentionally focused on relatively small systems with few components. For large systems, e.g., large power plants with multiple loops, a single PSM would not be sufficient to model the entire facility. Additional methods are needed to train and coordinate multiple PSMs representing separate parts of a plant. A potential approach may be to adopt a graphical architecture where PSMs represent transport between nodes.

Figure 13: A comparison of the spatial distribution of nominal and degraded residuals, \(r\) and \(r_{\mathrm{deg}}\), respectively, in Alg. 3. The residual values have been scaled by their maximum and minimum value.
## CRediT authorship contribution statement

**Akshay J. Dave**: Conceptualization, Methodology, Writing - Original Draft, Funding acquisition. **Richard B. Vilim**: Supervision, Methodology, Writing - Review & Editing.

## Acknowledgements

Argonne National Laboratory supports this work under the Laboratory Directed Research & Development Program, Project 2022-0077.
