Learning characteristic parameters and dynamics of centrifugal pumps
  under multiphase flow using physics-informed neural networks
Electrical submersible pumps (ESPs) are prevalently utilized as artificial
lift systems in the oil and gas industry. These pumps frequently encounter
multiphase flows comprising a complex mixture of hydrocarbons, water, and
sediments. Such mixtures lead to the formation of emulsions, characterized by
an effective viscosity distinct from that of the individual phases. Traditional
multiphase flow meters, employed to assess these conditions, are burdened by
high operational costs and susceptibility to degradation. To this end, this
study introduces a physics-informed neural network (PINN) model designed to
indirectly estimate the fluid properties, dynamic states, and crucial
parameters of an ESP system. A comprehensive structural and practical
identifiability analysis was performed to delineate the subset of parameters
that can be reliably estimated through the use of intake and discharge pressure
measurements from the pump. The efficacy of the PINN model was validated by
estimating the unknown states and parameters using these pressure measurements
as input data. Furthermore, the performance of the PINN model was benchmarked
against the particle filter method utilizing both simulated and experimental
data across varying water content scenarios. The comparative analysis suggests
that the PINN model holds significant potential as a viable alternative to
conventional multiphase flow meters, offering a promising avenue for enhancing
operational efficiency and reducing costs in ESP applications.
Felipe de Castro Teixeira Carvalho, Kamaljyoti Nath, Alberto Luiz Serpa, George Em Karniadakis
Learning characteristic parameters and dynamics of centrifugal pumps under multi-phase flow using physics-informed neural networks

###### Abstract

Electrical submersible pumps (ESP) are the second most used artificial lifting equipment in the oil and gas industry due to their high flow rates and boost pressures. They often have to handle multiphase flows, which usually contain a mixture of hydrocarbons, water, and/or sediments. Given these circumstances, emulsions are commonly formed. It is a liquid-liquid flow composed of two immiscible fluids whose effective viscosity and density differ from the single phase separately. In this context, accurate modeling of ESP systems is crucial for optimizing oil production and implementing control strategies. However, real-time and direct measurement of fluid and system characteristics is often impractical due to time constraints and economy. Hence, indirect methods are generally considered to estimate the system parameters. In this paper, we formulate a machine learning model based on Physics-Informed Neural Networks (PINNs) to estimate crucial system parameters. In order to study the efficacy of the proposed PINN model, we conduct computational studies using not only simulated but also experimental data for different water-oil ratios. We evaluate the state variable's dynamics and unknown parameters for various combinations when only intake and discharge pressure measurements are available. We also study structural and practical identifiability analyses based on commonly available pressure measurements. The PINN model could reduce the requirement of expensive field laboratory tests used to estimate fluid properties.

Electrical submersible pump Physics-informed neural networks Parameters estimation Identifiability analysis Multiphase flow Digital twin

## 1 Introduction

Multistage centrifugal pumps play a critical role in various industries, including oil and gas, water supply, chemical processes, and power generation. In the oil industry, these pumps often handle multiphase flows that contain a mixture of hydrocarbons, water, and/or sediment. Among the different multistage centrifugal pumps, since the introduction of Electrical Submersible Pumps (ESP) in the oil and gas industry in 1927, they have been the second most extensively used artificial lifting equipment [1]. The ESP is a multistage centrifugal pump driven by a submerged electrical motor. They are commonly used where a relatively high flow rate with high boosting pressure is necessary [2, 3, 4] e.g., in the petroleum industry. It is common in the petroleum industry to have two-phased liquid-liquid flows, which often
involve transporting oil and water together [5; 6]. In liquid-liquid (e.g., oil-water) flows, colloidal dispersions such as emulsions are commonly formed due to the chemical characteristics of the liquids. Emulsions occur when one of the liquids forms droplets (dispersed phase) and is immersed in the other liquid (continuous phase). Emulsions exhibit a non-Newtonian behavior that can cause instability in ESP operation [7].

Emulsions tend to separate the liquid phases due to oil-water or water-oil immiscibility. The separation time may vary from minutes to hours and may extend to days, months, or even years, depending on the stability of the emulsion. The droplet size distribution, chemical properties of the liquid phases, and the presence of surfactants determine the stability and type of emulsion (water-in-oil or oil-in-water) [8]. Furthermore, there may be a transition from one emulsion type to another, commonly called phase inversion [9; 10; 11]. Emulsion in the ESP system may be caused by turbulence and shear stress caused by the ESP, pipe system, and valve [12].

Furthermore, the effective viscosity of emulsion may considerably exceed the oil or water single-phase viscosity. It is affected by the oil and water phase individual viscosity, the amount of water in the oil-water mixture, the temperature, the droplet size, and distribution [13]. Therefore, monitoring fluid properties such as effective viscosity and density becomes crucial. These parameters affect the stability of the system and are important for developing control strategies and optimization of oil production. However, direct measurement of the fluid properties and system parameters next to the ESP becomes time-consuming and expensive. In this study, we considered a machine learning approach to identify these parameters using pressure measurements.

The problem of estimating unknown properties based on indirect measurements is referred to as the inverse problem and is common in various fields [14]. Additionally, multiple parameters may potentially provide satisfactory fits to a given data [15; 16; 17]. In this context, we conduct identifiability analysis to determine if a specific system parameter (or a set of parameters) can be accurately estimated uniquely (globally or locally) based on the available input and output states. The structural identifiability analysis assumes an ideal scenario, where the observed states are noise-free and the model itself is error-free [18; 17; 19]. This analysis, also known as prior identifiability, can be performed without actual experimental data.

Various methods have been developed to analyze a system's local and global structural identifiability. Chis _et al._[20] and Raue _et al._[21] conducted comparative studies of these methods found in the literature, focusing on biological systems. Castro and Boer [22] introduced a simple scaling method, which exploits the invariance of the equations under parameter scaling transformations. More recently, Dong _et al._[23] proposed a method based on differential elimination. It is important to note that structural identifiability alone does not account for the quantity and quality of available data or the numerical optimization algorithm employed [24; 19].

Several methods have been developed to estimate unknown parameters of ordinary differential equations (ODEs). The Nonlinear Least Squares (NLS) method is widely adopted due to its versatility and applicability to various ODE systems. However, it is computationally intensive, and the inaccuracy in numerical approximations of derivatives may pose challenges, particularly for stiff systems [25; 26; 27]. In addition to NLS, other methods, such as collocation methods [25], Gaussian process-based approaches [28], Bayesian methods [29], and the recently proposed two-stage approach using Neural ODEs [27], are a few alternative techniques for parameter estimation. The physics-informed neural network is one of the recent neural network-based methods for parameter estimation for ODEs/PDEs.

Raissi _et al._[30] proposed Physics-informed neural networks (PINNs), a new form of training of neural networks known, which is able to solve both forward and data-driven inverse problems. The method takes into account the physics of the problem while formulating the loss function. The physics of the problem is introduced in the loss function by way of minimizing the residue of the differential equations of the problem at some collocation points. Generally, the derivatives of the predicted unknown are calculated using automatic differentiation [31]. PINNs was further developed to different variants like conservative PINNs (cPINNs) [32], extended PINN (XPINN) [33], hp-VPINNs [34], Parareal PINN (PPINN) [35], Separable PINN [36] etc. Furthermore, McClenny and Braga-Neto [37] proposed a self-adaptive weights technique that automatically tuned the weights for different loss functions in a multi-loss function problem in PINN problem. In the present study, we consider self-adaptive weights for calculating the total loss. PINN and its variants were considered in inverse problems like unsaturated groundwater flow [38], diesel engine [39], supersonic flows [40] to name just a few.

Accurate modelling of the ESP system that considers the complexities of liquid-liquid two-phase flow presents significant challenges. To enhance tractability, a common approach involves employing a lumped-element model for the flow within the pump and pipes. However, this lumped model still
involves several unknown parameters, such as the fluid bulk modulus, effective viscosity of the emulsion, and varying equivalent resistance of the pipeline over the well's operational lifespan. Obtaining these parameters under the conditions prevalent in oil field extraction presents a significant challenge. Consequently, there is a need for techniques that can effectively obtain these parameters using limited and noisy data. Developing such techniques becomes crucial as they can optimize oil extraction and facilitate the implementation of control strategies within the ESP pumping system.

In this paper, we explore the possibility of using PINNs to estimate the important characteristic parameters (flow, pump, etc.) and predict the dynamics of state variables based on a limited set of known field state variables. We develop a PINN model for this problem and conduct computational studies with simulated and experimental data. First, we study the identifiability problem using structural and practical identifiability analyses. This helps us to determine whether a particular parameter and a combination of parameters are locally or globally identifiable. Once we obtain the possible set of unknown parameters, we formulate PINNs for the inverse problem. As discussed earlier, structural identifiability alone does not guarantee accurate identification of unknown parameters for noisy and experimental data. We conduct detailed numerical studies using simulated and experimental data with different water-to-oil ratios. The PINN solution could reduce the requirement for expensive field laboratory tests to obtain the flow parameters.

The remainder of this research paper is organized as follows: In Section 2, we discuss the methodology, we provide an overview of the experimental setup (SS2.1), the dynamic model of the ESP system (SS2.2), the structure of the PINN model (SS2.3), and the details regarding data acquisition and generation (SS2.4). We also discuss the structural identifiability analysis within this section (SS2.5). In Section 3, we present the results obtained from the structural local and practical identifiability analyses. Subsequently, in Section 4, we present the results of the PINN model for state and unknown parameter estimation. Finally, in Section 5, we present the conclusions drawn from this study, summarizing the key findings and discussing the results.

## 2 Methodology

This work consists of two parts: experimental and computational. The experimental part (Section 2.1) explains the materials and methods used to acquire data from a laboratory-scale ESP system. The computational part includes the definition of the PINN model in Section 2.3. We consider simulated and experimental data for the inverse problem. The data generation for the simulated study is discussed in Section 2.4. The identifiability analysis is discussed in Section 2.5.

### Experimental study

#### 2.1.1 Experimental setup

The experimental setup was designed to investigate the ESP performance under various flow conditions, including single-phase and water/oil two-phase flows. The main equipment is an eight-stage ESP model P100L 538 series from Baker Hughes(r), located downstream of the oil and water mixture point. A schematic diagram of the experimental setup is shown in Fig. 1, indicating all the components and dimensions of the piping system.

The ESP impellers have mixed flow geometry with an outer diameter of \(108\,\mathrm{mm}\) and have \(7\) blades. Pressure gauges and temperature transmitters were installed at the inlet and outlet of the ESP. They measure the suction and discharge pressure and temperature, respectively. The pressure transducers are capacitive and manufactured by Emerson Rosemount(tm), series 2088. The temperature transducers are four-wire resistance temperature detectors type PT100 with 1/10 DIN accuracy. The pump shaft torque and rotational speed are measured using the sensor model T21WN manufactured by HBM(r). Each pump on the experimental setup is driven by an electric three-phase induction motor supplied by a Variable Speed Driver (VSD) on each motor.

The fluids used in the experiments are water and a blend of mineral oil. The oil viscosity was characterized by a rotational rheometer model HAAKE MARS III (see Appendix A.1 for the viscosity measurements in different temperatures). The water fraction inside the oil phase is measured by the water cut meter model Nemko 05 ATEX 112 manufactured by Roxar(r). The Coriolis meters used to measure the density and the mass flow rates of the oil and water flow lines are model F300S355 manufactured by Micro Motion(r).

As the system operates in a closed loop, it tends to heat up and change the viscosity. Thus, the heat exchanger controls the fluid's temperature and the viscosity inside the flow line. The temperature is regulated by monitoring the ESP inlet temperature (\(T_{1}\)). This temperature measurement is also used
for calculating the working fluid viscosity. The heat exchanger (HE) has an independent water flow line comprised of a chiller, a heater, and a water tank, allowing it to cool and heat the fluid. Also, as the oil viscosity strongly depends on the temperature, it is possible to conduct tests with different viscosities with the same oil by adjusting the fluid's temperature.

The operational ranges and uncertainties of the instruments are presented in Table A.14. The specifications of the pumps, motors, variable speed drive, valve, tank, and heat control system are presented in Table A.13.

#### 2.1.2 Experimental procedure

In this section, we detail the experimental procedure utilized to evaluate the dynamic behavior of the ESP system when varying the ESP shaft angular velocity while operating under a two-phase flow of oil and water. For simplicity, when referring to ESP angular velocity, it refers to the ESP shaft angular velocity. The procedure is divided into two stages; the first stage is the start-up stage, and the objective of this stage is to obtain a stable water cut. The second stage is the transient acquisition stage, in which we collect the experiment data.

**Stage 1:**: Start-up:

1. Begin by slowly starting the oil pump and water pump at a low speed (e.g., \(62\,\mathrm{rad}\,\mathrm{s}^{-1}\)) to ensure system stability and prevent sudden pressure surges.
2. Gradually increase the rotation speed of the ESP while monitoring the suction pressure. Maintain the ESP angular velocity below \(105\,\mathrm{rad}\,\mathrm{s}^{-1}\) and ensure the suction pressure remains above \(100\,\mathrm{kPa}\) to prevent cavitation in the ESP first stages.

Figure 1: **Schematic diagram of the experimental setup:** A schematic diagram of the experimental setup indicates all the piping system’s components and dimensions. It has four distinct flow lines, an oil flow line, a water flow line, a two-phase flow line carrying the oil-water mixture, and a closed loop of water used in the heat exchanger (HE). The separation tank separates the oil and water phases from the emulsion by gravity and stores the oil and water phases. The fluids are pumped independently from the separator tank and operate in a closed loop. The twin-screw pump pumps the oil phase from the separator tank. Then, the oil phase flows through a shell and tube heat exchanger, a water cut meter (measures the water fraction inside the oil phase), and a Coriolis meter (measures the mass flow rate and density). A single-stage centrifugal pump pumps the water phase from the separator tank that flows through a Coriolis meter and a remotely controlled valve. Then, the oil and water phases mix in a “T” joint, forming an emulsion. Then, the emulsion is pumped by the ESP, flows through a remotely controlled valve, and returns to the separation tank.
3. Check the measured water cut and adjust the water pump angular velocity if necessary. Repeat this step until the desired water cut is achieved.
4. Stop the water pump.
5. Transient acquisition: 1. Increase the ESP angular velocity to the desired initial value. 2. Increase the oil pump angular velocity while closing the ESP downstream valve to maintain the suction pressure within the range of \(100\,\mathrm{kPa}\) to \(600\,\mathrm{kPa}\). The upstream pressure of the ESP should not exceed \(600\,\mathrm{kPa}\) for a safe operation. 3. Gradually increase the ESP angular velocity to reach the desired final value. 4. Increase the oil pump angular velocity while closing the downstream valve to maintain the suction pressure within the range of \(100\,\mathrm{kPa}\) to \(600\,\mathrm{kPa}\). 5. Return the ESP angular velocity to the desired initial speed. 6. Allow sufficient time for the suction, discharge pressure, and volumetric flow rate to stabilize. This stabilization period allows capturing only the dynamics of the ESP angular velocity change. 7. Start data acquisition. 8. Wait \(5\,\mathrm{s}\). 1. Adjust the ESP angular velocity on the VSD to the final desired value. 2. Allow another stabilization period for the suction, discharge pressure, and volumetric flow rate. This period allows the capture of the complete dynamic system response to the rotation change. 3. End data acquisition. 1. Return to the **Stage 1**:c

We considered the system stable when the coefficient of variation of the last \(20\,\mathrm{s}\) of measurements of the suction and discharge pressure and oil volumetric flow rate were smaller than \(0.8\,\mathrm{\char 37}\). These criteria are the same as used and described by Figueiredo _et al._[41]. Additionally, for the temperature, we considered a tolerance of \(\pm 0.25\,\mathrm{\SIUnitSymbolCelsius
Shannon sampling theorem for the maximum frequency specified by the pressure transmitter manufacturer. Then, the signal was downsampled to the same frequency.

### Electrical submersible pump system dynamic model

In this section, we describe the dynamic behavior of the ESP system by using a set of Ordinary Differential Equations (ODEs), derived from a bond graph representation of the system (presented in Fig. B.1 in Appendix B). The model is based on the works of Tanaka _et al._[42] and Higo _et al._[43], which provide bond graph models for centrifugal pumps and pipelines, respectively. In this study, we focused only on the mechanical and hydraulic domains of the system for simplicity. Therefore, the electrical motor was not considered in this model. Instead, we treated the shaft torque as an input to the system, decoupling the mechanical and electrical domains.

The ESP system is described by the state variable,

\[\mathbf{X}=\{Q_{p},\omega,Q_{1},Q_{2},P_{1},P_{2}\}, \tag{1}\]

where \(Q_{p}\) is the pump volumetric flow rate, \(\omega\) is the shaft angular velocity, \(Q_{1}\) and \(Q_{2}\) are the volumetric flow rate in the upstream and downstream pipeline respectively, \(P_{1}\) and \(P_{2}\) are the intake and discharge pressure respectively. The governing equations describing the dynamics of the system depend on the input vector \(\gamma(t)\). These are of the form

\[\dot{\mathbf{X}}=F(\mathbf{X},\gamma(t)). \tag{2}\]

The ESP system is governed by the following set of ODEs:

\[\frac{dQ_{p}}{dt}=\frac{\left(P_{1}-P_{2}+k_{3}\,\mu\,Q_{p}\right) A_{p}}{\rho\,L_{p}}+\frac{A_{p}\left(k_{1p}\,\omega\,Q_{p}+k_{2p}\,\omega^{2}+ k_{4p}\,{Q_{p}}^{2}\right)}{L_{p}}, \tag{3a}\] \[\frac{d\omega}{dt}=\frac{\gamma(t)-k_{1s}\,\rho\,{Q_{p}}^{2}-k_{2 s}\,\rho\,\omega\,Q_{p}-k_{3s}\,\mu\,\omega-k_{4s}\,\omega-k_{5s}\,\omega^{2}}{I_ {s}},\] (3b) \[\frac{dQ_{1}}{dt}=\frac{\left(k_{bd}\omega_{t}-Q_{1}\right)k_{bl} \,\mu+P_{in}-P_{1}-f_{f}(Q_{1},\mu,L_{u},d_{u})\,{Q_{1}}^{2}\,A_{u}}{\rho\,L_{ u}}-\frac{k_{u}\,{Q_{1}}^{2}}{2\,L_{u}\,A_{u}},\] (3c) \[\frac{dQ_{2}}{dt}=\frac{\left(P_{2}-P_{out}-f_{f}(Q_{2},\mu,L_{d}, d_{d})\,{Q_{2}}^{2}\right)A_{d}}{\rho\,L_{d}}-\frac{k_{d}\,{Q_{2}}^{2}}{2\,L_{d} \,A_{d}}-\frac{{Q_{2}}^{2}\,A_{d}}{L_{d}\,C_{v}(a)^{2}\,\rho_{0}},\] (3d) \[\frac{dP_{1}}{dt}=\frac{\left(Q_{1}-Q_{p}\right)B}{A_{u}\,L_{u}},\] (3e) \[\frac{dP_{2}}{dt}=\frac{\left(Q_{p}-Q_{2}\right)B}{A_{d}\,L_{d}}, \tag{3f}\]

where the parameters are described in detail in Appendix B, the function \(C_{v}(a)\) represents the valve coefficient as a function of the valve aperture (\(a\)), and \(\rho_{0}\) is the water density at \(15\,^{\circ}\mathrm{C}\). For simplicity, as we are not varying the valve aperture, we considered a single equivalent resistance coefficient obtained experimentally.

Considering the high viscosity of the emulsion, we assumed a laminar flow in the emulsion and oil flow lines. Therefore, the friction function \(f_{f}(Q,\mu,L,d)\) is derived using the Darcy-Weisbach equation and the friction factor for laminar flow, resulting in the following equation:

\[f_{f}(Q,\mu,L,d)=\frac{128\,L\,\mu}{\pi\,Q\,d^{4}}. \tag{4}\]

The ESP system model takes a single input variable, \(\gamma(t)\), representing the torque applied to the ESP shaft. In the experimental setup, the intake pressure of the twin-screw pump \(P_{in}\) and the pressure at the end of the pipeline \(P_{out}\) refers to the same separation tank that is open to the atmosphere. Therefore, for the sake of simplicity, these pressures are assumed to be constant and are set to \(0\,\mathrm{Pa}\).

It is important to mention that in both experimental and real-world conditions, various factors, such as turbulence and shear generated by the ESP, valve, and other devices installed in the flow line, contribute to emulsion formation, leading to changes in emulsion viscosity. Thus, the viscosity of the emulsion varies as the fluids traverse the system due to these effects. Additionally, the emulsion formation is influenced by the emulsion stability, chemical properties of the fluid, and the presence or addition of emulsifiers and surfactants. Moreover, in practical scenarios, temperature gradients along the pipeline significantly impact viscosity. In such cases, the assumption of a single viscosity is no longer valid, and appropriate assumptions must be made to account for these complexities.
In our specific case, the ESP operates within a closed-loop system. Consequently, the same emulsion continuously flows through the ESP and valves while the fluid temperature is controlled through a heat exchanger. Given the tank volume, the operating volumetric flow rate, temperature control measures, and the transient experiment duration, we assumed that viscosity remains approximately constant.

Although the Brinkman [44] model is widely used in the oil industry for emulsion effective viscosity, numerous alternative models for emulsion viscosity have been proposed (e.g., Einstein [45], Taylor [46] and Pal and Rhodes [47]). Recently, Bulgarelli _et al._[12] introduced a model specifically for effective emulsion viscosity inside ESPs. However, this model is exclusively applicable to ESPs and does not hold for pipelines. Thus, in this study, we employed the Brinkman [44] model to determine the emulsion viscosity for the entire system, it is expressed as:

\[\mu=\mu_{c}\left(\frac{1}{1-\Omega}\right)^{2.5}\qquad 0<\Omega<1, \tag{5}\]

where \(\Omega\) represents the water cut and \(\mu_{c}\) is the viscosity of the continuous phase

### Physics-informed neural networks

The physics-informed neural network (PINN) comprises of a neural network (NN) to approximate the solution of a given differential equation and a second part that is considered to encode the corresponding differential equation in the optimization of the network parameters. The differential equation is encoded in the loss function along with data loss with the help of automatic differentiation [31]. In this study, we consider a fully connected feed-forward deep neural network (DNN) to approximate the state variables (i.e. \(Q_{p}\), \(\omega\), \(Q_{1}\), \(Q_{2}\), \(P_{1}\), \(P_{2}\)). The DNN architecture comprises one input layer, \(h\) hidden layers, and one linear output layer and can be represented as follows:

\[\mathbf{y}_{0} =t\] (6a) \[\mathbf{y}_{i} =\sigma\left(\mathbf{W}_{i}\mathbf{y}_{i-1}+\mathbf{b}_{i}\right) i\ \forall\ 1\leq i\leq h-1\] Hidden layers (6b) \[\mathbf{\hat{y}} =\mathbf{y}_{h} =\mathbf{W}_{h}\mathbf{y}_{h-1}+\mathbf{b}_{h-1}\] Output layer (6c)

where \(\mathbf{W}\in\mathbb{R}^{m\times n}\) and \(\mathbf{b}\in\mathbb{R}^{m}\) are the weights and biases of the neural network. The \(n\) and \(m\) are the sizes of the previous and current layers, respectively. The function \(\sigma(.)\) is the activation function, which in this study we consider as a tanh\((.)\) function. For the remainder of this work, the weights and biases are referred to collectively as parameters of the neural network \(\mathbf{\theta}=\{\mathbf{W},\mathbf{b}\}\).

The schematic diagram of the PINN considered in this study is presented in Fig. 2. In this work, a single DNN approximates the set state variables of the ODEs presented in Section 2.2. This approach has demonstrated satisfactory results while maintaining a lower computational cost than training separate DNNs for each state variable.

The loss consists of three different losses, i.e., the physics loss (\(\mathcal{L}^{ode}(\mathbf{\theta},\mathbf{\Lambda},\mathbf{\lambda}_{r})\)), the data loss (\(\mathcal{L}^{data}(\mathbf{\theta},\mathbf{\lambda}_{d})\)) and the initial condition loss (\(\mathcal{L}^{ic}(\mathbf{\theta},\mathbf{\lambda}_{ic})\)). The total loss is given as,

\[\mathcal{L}(\mathbf{\theta},\mathbf{\Lambda},\mathbf{\lambda}_{d},\mathbf{\lambda}_{r},\mathbf{ \lambda}_{ic})=\mathcal{L}^{ode}(\mathbf{\theta},\mathbf{\Lambda},\mathbf{\lambda}_{r})+ \mathcal{L}^{data}(\mathbf{\theta},\mathbf{\lambda}_{d})+\mathcal{L}^{ic}(\mathbf{\theta},\mathbf{\lambda}_{ic}). \tag{7}\]

The physics loss is given by the weighted sum (self-adaptive weights) of physics loss due to individual equation,

\[\mathcal{L}^{ode}(\mathbf{\theta},\mathbf{\Lambda},\mathbf{\lambda}_{r}) =\sum_{s\in\Phi}m(\lambda_{r}^{s})\mathcal{L}^{ode}_{s},\qquad \Phi=\{Q_{p},\ \omega,\ Q_{1},\ Q_{2},\ P_{1},\ P_{2}\}, \tag{8a}\] \[=\sum_{s\in\Phi}m(\lambda_{r}^{s})\left[\frac{1}{sN^{ode}}\sum_{i =1}^{{}^{\prime}N^{ode}}s^{r_{i}}{}^{2}\right],\] (8b) \[=\sum_{s\in\Phi}m(\lambda_{r}^{s})\left[\frac{1}{sN^{ode}}\sum_{i =1}^{{}^{\prime}N^{ode}}\left(\frac{d\hat{y}_{s}}{dt}\bigg{|}_{t_{i}}-f_{s} \left(\hat{y}(t_{i};\mathbf{\theta}),\gamma(t_{i});\mathbf{\Lambda}\right)\right)^{2} \right], \tag{8c}\]

where \(\mathcal{L}^{ode}_{s}\), \(s\in\Phi\), \(\Phi=\{Q_{p},\ \omega,\ Q_{1},\ Q_{2},\ P_{1},\ P_{2}\}\) are the physic-informed loss for the differential equation describing the dynamics of the state variables (Equations (3a) to (3f)). \({}^{s}N^{ode}\) is the number of residual points for each of the equations, \(\lambda_{r}^{s}\) are the self-adaptive weights associated with each physics loss, and \(m(.)\) is a mask function which is considered as softplus function in the present study. \({}^{s}r_{n}\) are the residual of the equation. We discuss the calculation of the residual in Appendix C. In the present
study, the residual points for all the equations are considered the same. Thus, all \({}^{s}N^{ode}\) has the same points.

Similarly, the data loss is given by the weighted sum (self-adaptive weights) of data loss due to \(P_{1}\) and \(P_{2}\).

\[\mathcal{L}^{data}(\mathbf{\theta},\mathbf{\lambda}_{d}) =\sum_{s\in\phi}m(\lambda_{d})\mathcal{L}^{data}_{s},\qquad\phi= \{P_{1},\ P_{2}\} \tag{9a}\] \[=\sum_{s\in\phi}m(\lambda_{d})\left[\frac{1}{sN^{data}}\sum_{i=1} ^{{}^{s}N^{data}}\left[y_{s}(t_{i})-\hat{y}_{s}(t_{i};\mathbf{\theta})\right]^{2} \right], \tag{9b}\]

where \({}^{s}N^{data}\) is the number of points at which the measured data are available for each of the known quantities, \(y_{s}(t_{i})\) is the known value of the \(s\) quantity at time \(t_{i}\) and the corresponding approximated value from the neural network is \(\hat{y}_{s}(t_{i};\mathbf{\theta})\). Similar to physics loss, \(\lambda_{d}^{s}\) is the self-adaptive weight associated with each data loss, and \(m(.)\) is a mask function taken as the softplus function in the present study.

Figure 2: **Schematic representation of the Physics-informed neural network for the ESP system:** The deep neural network, shown in the red dashed-dotted rectangle on the left, is considered to approximate the solution of the ODEs system described in Section 2.2 (Equations (3a) to (3f)). The input to the neural network is time, denoted by \(t\), and the output is the ESP system states, highlighted in the green dotted rectangle of the figure. Each differential equation of the ESP system has residues at certain collocation points that must be minimized. We indicate them in the black dotted region on the right. The time derivative of each state (DNN output) is computed using automatic differentiation. The total loss, denoted as \(\mathcal{L}(\mathbf{\theta})\), includes data loss, physics (ODE) loss, and initial condition loss. The data loss (\(\mathcal{L}^{data}(\hat{y}-y)\)) is the loss between the DNN output and measured data, the physics (ODE) loss is the loss of the residue of the equations at the collocation points, and the initial condition loss is a loss between the initial condition and the DNN output at \(t=0\). \(\mathbf{\lambda}_{d}\), \(\mathbf{\lambda}_{r}\), and \(\mathbf{\lambda}_{ic}\) are the weights to the data loss, physics loss, and initial condition loss, respectively, while calculating the total loss. These may be fixed or adaptive. \(\mathbf{\Lambda}\) are the unknown parameters in the ODE system. The activation function of the DNN is represented by \(\sigma\).
The initial condition loss

\[\mathcal{L}^{ic}(\mathbf{\theta},\mathbf{\lambda}_{ic}) =\sum_{s\in\Phi}m(\lambda_{ic}^{s})\mathcal{L}_{s}^{ic},\qquad \Phi=\{Q_{p},\ \omega,\ Q_{1},\ Q_{2},\ P_{1},\ P_{2}\} \tag{10a}\] \[=\sum_{s\in\Phi}m(\lambda_{ic}^{s})\left[(y_{s}(t_{0})-\hat{y}_{c }(t_{0};\mathbf{\theta}))^{2}\right], \tag{10b}\]

where \(y_{s}(t_{0})\) is the known value of the \(s\) quantity at time zero, and the corresponding approximated value from the neural network is \(\hat{y}_{s}(t_{0};\mathbf{\theta})\). \(m(.)\) is a mask function taken as the softplus function in the present study.

The neural network architecture considered in the present study is \([1,20,20,20,6]\), where the input layer consists of one neuron that takes time \(t\) as input and is scaled between \([-1,1]\). In the output layer, we have six neurons that approximate the six state variables discussed in Section 2.2 (Equation (3a) to Equation (3f)). The neural network outputs pass through a transformation given by:

\[f(x)=\frac{(x+1)(x_{\text{max}}-x_{\text{min}})}{2}+x_{\text{min}} \tag{11}\]

where \(f(x)\) represents the state variable in the physical domain, and \(x\) is the output from the neural network. For known states (i.e., \(P_{1}\) and \(P_{2}\)), \(x_{\text{min}}\) and \(x_{\text{max}}\) are the minimum and maximum measured values, respectively. However, for the unknown states, the \(x_{\text{min}}\) and \(x_{\text{max}}\) are estimated by solving a system of equations at two specific time points \(t_{1}\) and \(t_{2}\). The system of equations, derived from simplifications of Equations (3a) and (3b), is given by:

\[\begin{cases}k_{1p}\,\rho\,\omega(t_{i})\,Q_{p}(t_{i})+k_{2p}\,\rho\,\omega(t_ {i})^{2}+k_{4p}\,\rho Q_{p}(t_{i})^{2}-(P_{2}(t_{i})-P_{1}(t_{i}))=0,\\ k_{1s}\,\rho\,Q_{p}(t_{i})^{2}-k_{2s}\,\rho\,\omega(t_{i})\,Q_{p}(t_{i})+\gamma (t_{i})=0,\qquad t_{i}\in\{t_{1},t_{2}\}\end{cases} \tag{12}\]

where \(k_{1p}\), \(k_{2p}\), \(k_{4p}\), \(k_{1s}\), and \(k_{2s}\) are the ESP system parameters, \(\rho\) represents the fluid density, and \(t_{i}\) denotes the time instant for the state variables (\(\omega\), \(Q_{p}\), \(P_{1}\), \(P_{2}\)) and system input (\(\gamma\)). Furthermore, the time points \(t_{1}\) and \(t_{2}\) are determined based on the pressure difference between \(P_{2}(t)\) and \(P_{1}(t)\) as follows:

\[t_{1}=\arg\min_{t}\left(P_{2}(t)-P_{1}(t)\right),\quad t_{2}=\arg\max_{t}\left( P_{2}(t)-P_{1}(t)\right). \tag{13}\]

Solving the system of equations at the time point \(t_{1}\) provides the estimation for \(x_{\text{min}}\), and at the time point \(t_{2}\) for \(x_{\text{max}}\). We would also like to mention that in order to solve the above Equation (12), we need to know the parameters (\(k_{1p}\), \(k_{2p}\), \(k_{4p}\), \(k_{1s}\) and \(k_{2s}\)) of the equation, which are unknown in this case. Thus, we assumed some initial value in the parameters. In this study, while solving the Equation (12), we considered \(\rho=931.51\,\mathrm{kg\,m^{-3}}\), which represents a water cut of \(50\,\%\), as an initial guess, and for the remaining pump parameters we considered a value of \(+15\%\) of the true value for all parameters. Additionally, we considered the same values obtained for \(Q_{p}\) for the \(Q_{1}\) and \(Q_{2}\) states (pipeline flow rates).

It is also worth noting that the state variables are at different scale (magnitude). To address this, we employ the transformation, described by Equation (11), to make the neural network's output of similar order for different state variables. This transformation can be of different forms, depending on the specific requirements of the problem. In the present study, the bounds \(x_{\text{min}}\) and \(x_{\text{max}}\) vary across the different experimental investigations as they depend on \(P_{2}(t)\) and \(P_{1}(t)\).

The trainable parameters (\(\mathbf{\theta}\), \(\mathbf{\Lambda}\), \(\mathbf{\lambda}_{d}\), \(\mathbf{\lambda}_{r}\) and \(\mathbf{\lambda}_{ic}\)) are optimized using the Adam optimizer. The automatic differentiation for the derivatives of the outputs with respect to input is evaluated using the Python library JAX, while the Adam optimization is carried out using the Optax library. As discussed earlier, we have considered self-adaptive weight [37]. Thus, we considered different optimizers (three) for the neural network parameters (weights and biases), unknown parameters, and self-adaptive weights. We discuss the number of training epochs and learning rate scheduler in detail in Appendix D for each numerical example. Further, the unknown parameters transformations will be discussed in Section 4.

### Data generation

In this study, we considered laboratory data collected from the experimental setup described in Section 2.1. Before considering the proposed PINN model for experimental data, we test out the model with simulated with and without noise. The simulated data are generated by solving the ODEs system presented in Section 2.2. This section will present the steps for generating the simulated data. We define the parameters and input used in the simulations. We try to generate the simulated data as close as to
the experimental data making the parameters and the input of the system near equal to the experimental setup, and are discussed hereafter.

For the simulated cases, in this study, we employed the Julia package DifferentialEquations.jl [48] to numerically solve the ESP system set of ODEs presented in Section 2.2. The simulated data were generated with \(\Delta t=$0.0001\,\mathrm{s}$\).

The model input is the torque (\(\gamma(t)\)), which we directly measured during each experimental investigation. After post-processing the torque signal, as described in Section 2.1.4, we observed that the torque signal still contained relatively high-frequency components, which would be unrealistic as an input for the ESP system. Therefore, we applied a Butterworth low-pass filter with a cut-off frequency of \(2\,\mathrm{Hz}\) to the torque signal to obtain a more representative input signal.

For fluid viscosity, we assumed a single value for the entire system. We calculated the average water cut and temperature observed during the corresponding experimental investigation and estimated the emulsion viscosity using the Brinkman [44] model. As for fluid density, we directly measured it using a Coriolis meter. Then, we calculated the average density obtained from the measurements conducted during each experimental investigation.

In the tuning step of a dynamic model, the parameters are adjusted to fit real-world observations. However, despite the tuning step, we expect discrepancies between the model predictions and measured data. They arise due to uncertainties in the data and the model's representation of the actual system, which may not fully capture the system complexities. Additionally, in inverse problems, minor variations or errors in the fixed parameter can significantly affect the desired parameter estimation.

Therefore, to accurately evaluate the capability of PINNs in handling measurement uncertainties while excluding potential influences of missing physics or the tuning steps, we incorporated Gaussian noise into the simulated signals \(P_{1}\) and \(P_{2}\). The noise values used in these signals were determined based on the manufacturer's specifications for the pressure transmitters' uncertainties, as detailed in Table A.14.

### Structural identifiability analysis

In the Section 2.2, we introduced a set of ODE that describes the ESP system dynamics that has \(26\) parameters and \(6\) states variables. Before trying to estimate the unknown parameters of the system using the PINN, it is important to assess whether or not the parameters can be uniquely determined from a given set of data.

Consider a dynamical system in the following format

\[\dot{\mathbf{x}}(t) =\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t),\mathbf{\Theta}), \tag{14}\] \[\mathbf{y}(t) =\mathbf{h}(\mathbf{x}(t),\mathbf{u}(t),\mathbf{\Theta}), \tag{15}\]

where \(\mathbf{x}(t)\) is an \(m\)-dimensional state vector, \(\mathbf{u}(t)\) is an \(n\)-dimensional input signal, \(\mathbf{y}(t)\) is an \(r\)-dimensional output signal or the measurable output, and \(\mathbf{\Theta}\) is \(k\)-dimensional vector of parameters. A parameter set \(\mathbf{\Theta}\) is said to be structurally globally identifiable if the following property holds:

\[\mathbf{h}(\mathbf{x}(t),\mathbf{u}(t),\mathbf{\Theta})=\mathbf{h}(\mathbf{x}(t),\mathbf{u}(t),\mathbf{\beta })\implies\mathbf{\Theta}=\mathbf{\beta}, \tag{16}\]

where \(\mathbf{\beta}\) is a \(k\)-dimensional vector of parameters. Furthermore, if the property in Equation (16) holds within a neighborhood of \(\mathbf{\Theta}\), it is referred to as structurally locally identifiable.

Thus, the identifiability property serves as a prerequisite for practical parameter estimation. Additionally, as stated by Daneker _et al._[49], if a parameter is locally identifiable, it implies that the search range for that parameter should be limited before attempting its estimation. On the other hand, for globally identifiable parameters, there is no need to define a search range.

## 3 Results and discussions for identifiability analysis

### Local identifiability

We utilized the method of differential elimination for dynamical models via projections proposed by Dong _et al._[23] to evaluate the structural identifiability of the ESP system model. It is important to note that our analysis focused specifically on assessing local identifiability due to out-of-memory issues when attempting to test for structural global identifiability. Additionally, our investigation considered the suction and discharge pressure (\(P_{1}\) and \(P_{2}\)) as available measurements, as they are typically measured in ESP deployments within oil fields.

From Tables B.1 and B.2, we observe that the ESP system model initially contains a total of \(26\) parameters. However, it is worth noting that \(8\) of these parameters are geometrical parameters, such
as pipe diameter, cross-sectional area, and pipeline length. They can be readily obtained, and they were fixed for both the upstream and downstream pipelines. Furthermore, \(3\) parameters refer to the twin-screw pump, which serves as a pressure booster for the flow line. In actual oil field extraction, the pressure is a characteristic of the well. Thus, we reduce the number of unknown parameters to \(15\). Subsequently, we assessed the local identifiability of the ESP system model based on these \(15\) parameters. The corresponding results are displayed in the first row of Table 2.

The ESP system model with \(15\) unknown parameters is found to be structurally locally unidentifiable, and it is impossible to estimate all of the parameters simultaneously. As discussed by Daneker _et al._[49], there are several approaches to address the identifiability issue, such as fixing specific parameters or introducing additional measured variables. For simplicity, when we refer to fixing a parameter, we consider that the parameter is known. In this study, we opted to work with the available measurements from the actual field and consider only the suction and discharge pressures (\(P_{1}\) and \(P_{2}\)). Therefore, we decided to fix the shaft inertia (\(I_{s}\)). Despite this adjustment, the model remains structurally locally unidentifiable, with only \(k_{3s}\) and \(k_{4s}\) being locally unidentifiable. The corresponding results are shown in the second row of Table 2.

Finally, we solved the issue of the model being locally unidentifiable by fixing both \(k_{3s}\) and \(k_{4s}\). The results are shown in the third row of Table 2. We found that fixing only one of them did not solve the problem. In summary, the ESP system model can only be structurally locally identifiable with the suction and discharge pressure measurements if the shaft inertia (\(I_{s}\)), viscous damping coefficient (\(k_{4s}\)), and disk friction coefficient (\(k_{3s}\)) are known.

### Practical identifiability

In section Section 3.1, we discussed the circumstances under which the ESP system is structurally locally identifiable. However, it is worth noting that this analysis was conducted assuming that the measured variables have no noise and that the model is error-free. Hence, we must also assess if we can estimate parameters accurately based on the available data's quantity and quality. This analysis is known as practical or posterior identifiability analysis, and it is possible that a structurally identifiable system may not be practically identifiable.

The practical identifiability can be assessed using either Monte Carlo simulations or sensitivity analysis. The Monte Carlo approach rigorously determines the practical identifiability of the model but comes with a high computational cost due to the requirement of multiple model fits. On the other hand, sensitivity analysis offers a faster computation method and provides information about the correlation structure among the parameters. This correlation structure can guide the fixing of parameters when practical identifiability is not achieved [19; 50].

In this study, we considered only the sensitivity approach to practical identifiability. Then, we used the Fisher information matrix (FIM) to compute the correlation matrix of all parameters to determine their practical identifiability. We estimate the sensitivities of the ESP system model with respect to the parameters with the Julia package DiffEqSensitivity.jl. For the sensitivity analysis, we considered the simulated case, and similarly to Daneker _et al._[49], we considered a noise level of the measurements of \(1\,\%\).

Then, to obtain the correlation matrix, we first estimate the covariance matrix (\(\mathbf{C}\)) which can be approximately obtained from the FIM by

\[\mathbf{C}=FIM^{-1}. \tag{17}\]

We then estimate the correlation matrix from (\(\mathbf{C}\)) with

\[\begin{cases}r_{ij}=\frac{C_{ii}}{\sqrt{C_{ii}}\,C_{ij}},&\text{if }i\neq j,\\ r_{ij}=1,&\text{if }i=j.\end{cases} \tag{18}\]

\begin{table}
\begin{tabular}{c c c c c c c c c c c c c c} \hline \hline \(k_{1p}\) & \(k_{2p}\) & \(k_{3p}\) & \(k_{4p}\) & \(k_{1s}\) & \(k_{2s}\) & \(k_{3s}\) & \(k_{4s}\) & \(k_{5s}\) & \(I_{s}\) & \(\rho\) & \(B\) & \(k_{u}\) & \(k_{d}\) \\ \hline \(\times\) & \(\times\) & ✓ & ✓ & \(\times\) & \(\times\) & \(\times\) & \(\times\) & \(\times\) & \(\times\) & ✓ & ✓ & ✓ & ✓ & ✓ & ✓ \\ ✓ & ✓ & ✓ & ✓ & ✓ & ✓ & \(\times\) & \(\times\) & ✓ & \(-\) & ✓ & ✓ & ✓ & ✓ & ✓ \\ ✓ & ✓ & ✓ & ✓ & ✓ & ✓ & \(-\) & \(-\) & ✓ & \(-\) & ✓ & ✓ & ✓ & ✓ & ✓ \\ \hline \hline \end{tabular}
\end{table}
Table 2: **Local structural identifiability results of the ESP system model, considering that the intake pressure (\(P_{1}\)) and discharge pressure (\(P_{2}\)) are known.** The symbols \(\times\) indicate locally unidentifiable parameters, ✓ indicate locally identifiable parameters, and \(-\) represent fixed parameters.
When \(|r_{ij}|\) is close to \(1\), the parameters \(i\) and \(j\) are strongly correlated and cannot be individually estimated. Therefore, the parameters are practically unidentifiable [19, 50]. The absolute of the correlation matrix of the parameters structurally locally identifiable of the ESP system is shown in Fig. 2(a). For the practical identifiability analysis, we considered the simulated case for the first investigation of the test matrix, Table 1.

The correlation analysis shown in Fig. 2(a) clearly demonstrates a significant correlation among multiple parameters. Notably, the fluid density (\(\rho\)) and the pipeline resistances (\(k_{u}\) and \(k_{d}\)) exhibit a strong correlation. Consequently, although the initial analysis indicated structural local identifiability of these parameters, they are not practically identifiable within the context of the conducted investigation.

As many parameters are strongly correlated, we need to fix some parameters to make the system practically identifiable. We set that the main variables of interest are the fluid parameters and the pipeline resistances since the fluid parameters (\(B\), \(\mu\), and \(\rho\)) change within the water cut, which changes with the well's life, and the pipeline resistance (\(k_{u}\) and \(k_{d}\)) changes with the actuation of valves and wax deposition for instance. Thus, we preferably fixed the pump and shaft parameters. We started by fixing the parameter \(k_{1p}\), then the \(k_{3s}\), \(k_{2p}\) and \(k_{3p}\) successively. The resultant parameter correlation matrix is presented in Fig. 2(b).

It is noticeable from Fig. 2(b) that after fixing some parameters, they are less correlated. However, some parameters are still correlated, such as the viscosity (\(\mu\)) and the pump equivalent resistance (\(k_{4p}\)). However, this correlation is not as strong as the correlation between density (\(\rho\)) and pipeline resistance (\(k_{u}\)) observed in Fig. 2(a). In Fig. 2(b), the correlation coefficient between \(\mu\) and \(k_{4p}\) is \(0.9880\), while in Fig. 2(a), the correlation coefficient between \(\rho\) and \(k_{u}\) is \(0.9994\).

## 4 Results and discussions for inverse problem using PINN

From the structural identifiability analysis presented in Section 3, we have defined three cases to evaluate the PINN for parameter and state estimation in the ESP system. Consistent with Section 3, we considered only the suction and discharge pressures as the measured variables. The three cases are as follows:

1. Flow parameters: In this case, we assumed as known all parameters of the ESP system model except for the flow properties, namely, density (\(\rho\)), viscosity (\(\mu\)), and the bulk modulus (\(B\)).
2. Flow, pipeline, and impeller parameters: Based on the local identifiability analysis results presented in the first row of Table 2, we assumed as known the parameters marked as locally unidentifiable and estimate the remaining identifiable ones. The unknown parameters considered are the flow properties (\(B\), \(\mu\), \(\rho\)), the upstream pipeline equivalent resistances (\(k_{u}\) and \(k_{d}\)) and the pump parameters \(k_{3p}\) and \(k_{4p}\).

Figure 3: **Absolute values of the estimated correlation matrices for ESP system model parameters. The matrices are derived based on the FIM. The matrix indicates whether parameters are identifiable or not. Parameters with absolute correlation values near 1 are practically unidentifiable. The Fig. 2(a) represents twelve structurally locally identifiable parameters, while the Fig. 2(b) focuses on eight parameters, particularly keeping fluid properties and pipeline resistances unknown.**
* Flow, pipeline, impeller, and shaft parameters: Considering the findings from the practical identifiability analysis outlined in Section 3.2, we estimate the parameters depicted in Fig. 2(b) while assuming as known the remaining parameters. The unknown parameters are fluid properties (\(B\), \(\mu\), \(\rho\)), the upstream pipeline equivalent resistances (\(k_{u}\) and \(k_{d}\)), the pump parameter \(k_{4\text{p}}\), and the shaft parameters \(k_{1s}\) and \(k_{5s}\).

For each case, we will evaluate the performance and effectiveness of PINN in estimating the unknown parameters and non-measured system states in three distinct data scenarios:

1. Simulated data: This scenario represents an ideal condition without any disturbances or modeling errors, allowing us to assess the baseline performance of the proposed method.
2. Simulated data with Gaussian noise: In this scenario, we introduce Gaussian noise to the simulated data to evaluate the sensitivity of the proposed method to noisy data.
3. Experimental data: This scenario is closer to actual oil field conditions, incorporating noisy measurements, modeling errors, and missing information about certain parameters.

Furthermore, for each study, we trained the PINN \(30\) times to evaluate the impact of the neural network weights initialization. The results are presented and discussed in each case for the two experimental conditions defined in Table 1. Fig. 4 presents a schematic representation of the cases and scenarios evaluated.

During the training of the PINN, we observed that the significant difference in magnitude between the states, with pressure in the order of \(10^{5}\) and flow rates in the order of \(10^{-2}\), affected the training and tuning process despite the scaling in the physics loss and the usage of self-adaptive weights on the loss function. Therefore, we converted the units of the states from SI units to \(\mathrm{m}^{3}\,\mathrm{h}^{-1}\) for the volumetric flow rates and to \(\mathrm{MWC}\) (meters of water column) for the pressure. After obtaining the results, we converted them back into SI units.

Furthermore, the transformation of the unknown parameters contributes to the convergence and accuracy of the PINN by restricting the search space to the local neighborhood of the parameter. We first tried to keep the linear scaling for all parameters and cases during the development. However, in Case 3, we could not achieve satisfactory results without restricting the parameter search space using bounded transformations such as \(\mathrm{tanh}\). Also, we took advantage of that in our experiments. It is not possible to have a density higher than the water. Then, we set an upper bound for it. It helped to improve the results in all cases. The transformation scheme for each parameter is shown in Table 3, and the scaling values are presented in Table 4.

Figure 4: **Cases and scenarios evaluated.** As presented in Section 2.1.3, two experimental investigations (blue) are performed, considering different water cuts and initial angular velocities. Each experimental investigation is analyzed under three different data sources (scenarios): simulation data, simulation data with added Gaussian noise, and experimental data collected from the instruments (green). For each scenario, three sets of unknown parameters, denoted as cases, are evaluated (red). Thus, we have \(9\) different results for each experimental investigation.
### State estimation results

We begin by assessing the state estimation capabilities of the PINN when applied to the ESP system model. The PINN will estimate the remaining unknown system states using the suction and discharge pressure signals as known variables (\(P_{1}\) and \(P_{2}\)). In the simulated cases, the unknown states include the pipeline and ESP volumetric flow rates (\(Q_{1}\), \(Q_{2}\), \(Q_{p}\)) and the ESP angular velocity (\(\omega\)). However, in the experimental case, measurements of the downstream pipeline and ESP volumetric flow rates (\(Q_{2}\) and \(Q_{p}\)) are unavailable, and therefore, they were excluded from the accuracy analysis. We would also like to mention that these results are from inverse problems with unknown parameters discussed above, not forward problems. First, we present the state estimation results for all the cases considered. Then, we discuss the prediction of the unknown parameters for all the cases in Section 4.2.

It is worth mentioning that in PINN, it is common practice to use second-order optimization algorithms, such as L-BFGS, after training with the first-order algorithm. This additional step usually aims to fine-tune the PINN's performance. However, we observed that using L-BFGS after Adam adversely affected the parameter estimation in the noisy and experimental case. Therefore, we decided to rely solely on the first-order optimization algorithm, as our focus is on unknown parameter estimation in more realistic scenarios. It is important to note that the simulated cases serve as a baseline performance for assessing the PINN's performance, and we expect that performance will be comparatively reduced under the more challenging conditions posed by experimental data and noise.

#### 4.1.1 State estimation results for simulated data

As a baseline performance assessment, we begin with the simulated scenario. The predicted states from the neural network are shown in Fig. 5, while the Mean Absolute Percentage Error (MAPE) corresponding to the prediction of different states is shown in Table 5. The MAPE is defined as

\[\text{MAPE}=\frac{1}{N^{map}}\sum_{i=1}^{N^{map}}\left|\frac{y_{s}(t_{i})-\bar {y}_{s}(t_{i};\mathbf{\theta})}{y_{s}(t_{i})}\right|, \tag{19}\]

where \(N^{map}\) is the total number of temporal data points. We considered the same sample rate as the experimental data described in Section 2.1.4, \(10\,\mathrm{Hz}\) (\(\Delta t=0.1\) sec). The variable \(y_{s}(t_{i})\) represents the true value of the state \(s\) at time \(t_{i}\), and \(\bar{y}_{s}(t_{i};\mathbf{\theta})\) denotes the mean of the predicted value of the state \(s\) from the PINN across 30 different realization of the results. These results are evaluated with different initialization of the parameters of the network and unknown parameters. For the sake of brevity, we

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline
**Case** & \(B\) & \(\mu\) & \(\rho\) & \(k_{u}\) & \(k_{d}\) & \(k_{3p}\) & \(k_{4p}\) \\ \hline
**Case 1** & \(1\times 10^{9}\) & \(1\) & \(1000\) & & & & \\
**Case 2** & \(1\times 10^{9}\) & \(1\) & \(1000\) & \(3\times 10^{2}\) & \(1\times 10^{1}\) & \(-1.0\times 10^{7}\) & \(-1.0\times 10^{2}\) \\
**Case 3** & & & \(1000\) & & & & \\ \hline \hline \end{tabular}
\end{table}
Table 4: **Scaling constants for the unknown parameters. The table presents the scaling constants used in linear parameter transformations and the scaling parameter for the liquid bulk modulus (\(B\)) and density (\(\rho\)) transformations.**

\begin{table}
\begin{tabular}{l c c c c c c c c} \hline \hline
**Case** & \(B\) & \(\mu\) & \(\rho\) & \(k_{u}\) & \(k_{d}\) & \(k_{3p}\) & \(k_{4p}\) & \(k_{1s}\) & \(k_{5s}\) \\ \hline
**Case 1** & \((S_{p}(x)+0.9)\Lambda\) & \(L(x)\) & \(S_{m}(x)\Lambda\) & & & & & \\
**Case 2** & \((S_{p}(x)+0.9)\Lambda\) & \(L(x)\) & \(S_{m}(x)\Lambda\) & \(L(x)\) & \(L(x)\) & \(L(x)\) & & \\
**Case 3** & \(\Lambda_{\text{sc}}(x)\) & \(\Lambda_{\text{sc}}(x)\) & \(S_{m}(x)\Lambda\) & \(\Lambda_{\text{sc}}(x)\) & \(\Lambda_{\text{sc}}(x)\) & & \(\Lambda_{\text{sc}}(x)\) & \(\Lambda_{\text{sc}}(x)\) & \(\Lambda_{\text{sc}}(x)\) \\ \hline \hline \end{tabular}
\end{table}
Table 3: **Transformations of the unknown parameters for each case evaluated. The transformations applied to the unknown parameters are detailed. The softplus function \(S_{p}(\cdot)\) and the \(S_{m}(\cdot)\), defined as \(S_{m}(x)=x-\text{softplus}(x)+1\), are utilized along with a scaling constant \(\Lambda\). The \(L(\cdot)\) is a linear transformation given by \(L(x)=x\Lambda\), which considers only a scaling constant \(\Lambda\). For Case 3, the transformation for an unknown parameter is expressed as \(\Lambda_{\text{sc}}(x)=\Lambda_{\text{true}}\left(\tanh(x)\alpha+1\right)\), where \(x\) represents the PINN-estimated value of the unknown parameter, \(\Lambda_{\text{true}}\) is the parameter true value, and \(\alpha\) denotes the span percentage around the true value. For most parameters, we set the span to \(\pm 50\%\) (\(\alpha=0.5\)), except for the \(B\) parameter, which spans \(\pm 15\%\) (\(\alpha=0.15\)).**
restrict our discussion to the state prediction results and the MAPE table for the first experimental investigation. Further results from the second investigation are presented in Appendix E. We have shown the mean value of the predicted parameters and standard deviation for all the cases in Tables 11 and 12.

The comparison between the predicted dynamics of the states by the PINN model and the actual values, as shown in Fig. 5, demonstrates a good overall agreement across all states and cases. Notably, a small error is observed for the \(P_{1}\) and \(P_{2}\) states in all cases. However, a slightly larger error is observed for the volumetric flow rate states (\(Q_{1}\), \(Q_{2}\), and \(Q_{p}\)) and the ESP angular velocity state (\(\omega\)). Among the different cases, Case 1 (red line) exhibits the best agreement across the entire simulation range in both investigations. This is also observed in Table 5, where the Case 1 presented the lower MAPE in all states except \(P_{1}\).

In the first investigation, the flow rates in Case 2 (black line) and Case 3 (green line) show a reasonable agreement with the simulation approximately from the \(6^{th}\) second. However, a slight offset

\begin{table}
\begin{tabular}{l l l l l l l} \hline  & \(P_{1}\) & \(P_{2}\) & \(Q_{1}\) & \(Q_{2}\) & \(Q_{p}\) & \(\omega\) \\ \hline Case 1 & \(0.142\,\%\) & \(0.008\,\%\) & \(0.004\,\%\) & \(0.006\,\%\) & \(0.007\,\%\) & \(0.025\,\%\) \\ Case 2 & \(0.147\,\%\) & \(0.009\,\%\) & \(0.039\,\%\) & \(0.038\,\%\) & \(0.039\,\%\) & \(0.872\,\%\) \\ Case 3 & \(0.932\,\%\) & \(0.010\,\%\) & \(0.041\,\%\) & \(0.040\,\%\) & \(0.041\,\%\) & \(0.618\,\%\) \\ \hline \end{tabular}
\end{table}
Table 5: **MAPE for the simulated scenario.** The table presents the results for the simulation data considering three different cases of unknown parameters. The MAPE values represent the error of the simulation results for the states: \(P_{1}\), \(P_{2}\), \(Q_{1}\), \(Q_{2}\), \(Q_{p}\), and \(\omega\). As the number of unknown parameters increased, from Case 1 to Case 3, the MAPE values also increased, indicating a slight performance loss. However, the MAPE values remain relatively small, and visual inspection of Fig. 5 reveals a reasonable agreement with the actual values.

Figure 5: **Comparison of predicted and true states for simulated scenario of the first experimental investigation with unknown parameters.** The blue line indicates the true values, while the \(\times\) markers denote the data points used for training the PINN. The red, black, and green lines represent the mean values of the predicted states calculated at each time instant across the \(30\) trained PINNs for Cases 1 to 3, as defined earlier in this section (Section 4). Across all cases, the training dataset consists of \(30\) data points for the data loss in \(P_{1}\) and \(P_{2}\) and \(100\) collocation points for the physics loss.
from the actual value is observed from the beginning to the \(6^{th}\) second. As for the second investigation (Fig. E.1), a slightly higher offset is observed in the flow rate states and in the ESP angular velocity across the entire simulation time for Case 3, while for Case 2, we no longer observe the offset.

#### 4.1.2 State estimation results for simulated data with noise

To further evaluate the robustness of the PINN model, we extend our analysis to include a simulated scenario with Gaussian noise. This section aims to assess the performance of the PINN model when subjected to noise alone without considering missing physics or errors in the model parameter estimation. For each PINN realization, a different noise realization is also considered. Similarly to the previous section, we will focus on presenting the result corresponding to the first experimental investigation, while the result for the second experimental investigation can be found in Appendix E for reference. The MAPE for the different predicted states are shown in Table 6 and the output of the predicted states from the neural network along with the addition of noise on the pressure signals are shown in Fig. 6.

The comparison between the predicted dynamics of the states by the PINN and the actual values, as depicted in Fig. 6, demonstrates a good overall agreement all states and investigations for almost all cases. However, in Case 2 (black line), the PINN was unable to accurately estimate the flow rate states (\(Q_{1}\), \(Q_{2}\) and \(Q_{p}\),). They exhibit relatively higher errors, however, the model is able to predict the general shape of the time series. Remarkably, these higher errors were observed exclusively in the first experimental investigation (Fig. 7) (higher water cut), whereas in the second investigation (Fig. E.3) (with low water cut), Case 2 demonstrated performance similar to that of Case 1 (red line).

Notably, minor errors are observed in the \(P_{1}\) and \(P_{2}\) states for all cases and investigations. Specifically, in the discharge pressure state (Fig. 6b), we can observe that the neural network did not overfit and successfully captured the signal trend and magnitude. Similarly to the simulated data without noise (Section 4.1.1), Case 1 (red line) exhibits the best agreement across the entire simulation range in both investigations. Additionally, as shown on Table 6, Case 1 had a lower MAPE than the other cases for the unknown states and a higher MAPE for the known state \(P_{2}\).

Despite the addition of noise, the state estimation for the Case 1 and 3 provides good accuracy without significant errors. These results were similar to those obtained in noise-free conditions, with a slight increase in MAPE values. In Case 3, there is a noticeable offset in the volumetric flow rates (\(Q_{1}\), \(Q_{2}\) and \(Q_{p}\)) approximately after the \(8^{\text{th}}\) second, which is different from the observations made in Section 4.1.1, where the offset occurred only in the beginning. Additionally, the angular velocity, in this case, remains closely aligned with the actual values throughout the simulation, exhibiting no offset. It is observed that results for Case 2 deviate from the true values, particularly, for flow rates (\(Q_{1}\), \(Q_{2}\), \(Q_{p}\)). It may be possible that, though the Case 2 is structurally identifiable, it is not practical identifiable. However, it is observed that the angular velocity of the ESP (\(\omega\)) was not significantly affected as that of the flow rates.

\begin{table}
\begin{tabular}{l l l l l l l} \hline \hline  & \(P_{1}\) & \(P_{2}\) & \(Q_{1}\) & \(Q_{2}\) & \(Q_{p}\) & \(\omega\) \\ \hline Case 1 & \(0.180\,\%\) & \(0.073\,\%\) & \(0.017\,\%\) & \(0.017\,\%\) & \(0.017\,\%\) & \(0.027\,\%\) \\ Case 2 & \(0.165\,\%\) & \(0.021\,\%\) & \(0.455\,\%\) & \(0.456\,\%\) & \(0.455\,\%\) & \(2.815\,\%\) \\ Case 3 & \(0.182\,\%\) & \(0.034\,\%\) & \(0.048\,\%\) & \(0.045\,\%\) & \(0.048\,\%\) & \(0.159\,\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 6: **MAPE for the simulated scenario with added noise.** The table presents the MAPE for states prediction in the simulated case with added noise data, considering three cases of unknown parameters. Similarly to Table 5, increasing the number of unknown parameters (Case 1 to Case 3) leads to higher MAPE values for all states, indicating a performance loss. Nevertheless, the MAPE values remain relatively small, and visual inspection of Fig. 6 shows a reasonable agreement with the actual values, except for Case 2.
#### 4.1.3 State estimation results for experimental data

In this section, we assess the performance of the PINN when subjected to experimental data, which is a more representative scenario of real-world oil field conditions. In this case, noise in data and model uncertainties exist, posing a more challenging scenario for the PINN model. Similarly to the previous sections, the results from the first experimental investigation are provided in this section, while results from the second investigation can be found in Appendix E. The MAPE values for these predictions are presented in Table 7. Meanwhile, the predicted outputs of the neural network are shown in Fig. 7.

By comparing the PINN-predicted dynamics with the actual values (Fig. 7), a good overall agreement is evident across all states and cases, even in the presence of noise and model errors. However, similarly to the noisy data scenario, in Case 2 (black line), the PINN was unable to accurately estimate the flow rate (\(Q_{1}\)) and angular velocity (\(\omega\)) states. While the predicted states follow the signal's shape, they exhibit relatively high errors. Remarkably, these higher errors were observed exclusively in the first experimental investigation (Fig. 7) (high water cut), whereas in the second investigation (Fig. E.3) (low water cut), Case 2 demonstrated performance similar to that of Case 3 (green line). These high errors are observed in Table 7, where the Case 2 values for \(Q_{1}\) and \(\omega\) had the highest MAPE values.

Notably, for the \(P_{1}\) state (Fig. 6(a)), a minor deviation in the PINN prediction's shape is noticeable when compared to the experimental data. This difference could be attributed to the influence of the

Figure 6: **Comparison of predicted and true states for the simulated case with noise of the first experimental investigation and unknown parameters. In Fig. 5(a) and Fig. 5(b), the blue line indicates the signal with added Gaussian noise. In other figures, the same blue line represents the noise-free simulated state. This differentiation aids in evaluating the PINN’s accuracy by comparing its predictions with the actual state values and also present the signal used for PINN training. Each PINN training used a different random seed. Therefore, the Gaussian noise in the signal varied for every PINN training and case. Specifically, the noisy signals for \(P_{1}\) and \(P_{2}\) in Fig. 5(a) and Fig. 5(b) correspond to the first training for Case 1, and they were used as a reference. The added Gaussian noise considered the uncertainty from the pressure transducers as specified in Table A.14. For all cases, the training dataset consisted of \(37\) data points for data loss in \(P_{1}\) and \(P_{2}\) and \(100\) collocation points for physics loss. The plots show the mean value of results for 30 realization of PINNs. The noisy data shown in the plots are one of noisy data shown as representative sample.**
fluid's bulk modulus, which affects the system's settling time. In the parameter estimation section, we will discuss the challenges encountered in estimating the bulk modulus for the experimental data, as it was not experimentally measured. Furthermore, in the discharge pressure state (Fig. 6(b)), it is clear that the neural network did not overfit and successfully captured the signal trend and magnitude, including the signal overshoot between \(6\,\mathrm{s}\) to \(12\,\mathrm{s}\), which was considerably smaller in the simulations.

Additionally, a systematic error is noticeable at the beginning of the \(\omega\) state (Fig. 6(d)) for Cases 1 and 3. Despite these deviations, the predicted \(\omega\) signal aligns well with the experimental data regarding shape and magnitude. When comparing the simulated and measured angular velocity (\(\omega\)), we also observed a systematic error, indicating a difference between the actual and the obtained model parameters. However, the PINN could reasonably estimate the angular velocity state despite this differ

\begin{table}
\begin{tabular}{l l l l l} \hline  & \(P_{1}\) & \(P_{2}\) & \(Q_{1}\) & \(\omega\) \\ \hline Case 1 & \(4.517\,\mathrm{\char 37}\) & \(0.245\,\mathrm{\char 37}\) & \(0.137\,\mathrm{\char 37}\) & \(1.065\,\mathrm{\char 37}\) \\ Case 2 & \(3.828\,\mathrm{\char 37}\) & \(0.243\,\mathrm{\char 37}\) & \(1.560\,\mathrm{\char 37}\) & \(8.290\,\mathrm{\char 37}\) \\ Case 3 & \(2.963\,\mathrm{\char 37}\) & \(0.252\,\mathrm{\char 37}\) & \(0.218\,\mathrm{\char 37}\) & \(1.300\,\mathrm{\char 37}\) \\ \hline \end{tabular}
\end{table}
Table 7: **MAPE for the experimental data scenario.** The table presents the MAPE values for state predictions using experimental data across the three cases with varying unknown parameters. Among the cases evaluated, Case 1 shows the best agreement regarding the upstream volumetric flow rate (\(Q_{1}\)) and the ESP angular velocity \(\omega\). However, similarly to the noisy data scenario, Case 2 presented errors for \(\omega\) and \(Q_{1}\) that are considerably higher than the other states, as evident in Figs. 6(c) and 6(d). Nonetheless, the MAPE values remain relatively small for the other cases and states, and visual inspection of Fig. 6 shows a reasonable agreement with the actual values.

Figure 7: **Comparison of predicted and true states for the experimental case of the first experimental investigation with unknown parameters.** Due to constraints in our experimental setup, we could only directly measure the states \(Q_{1}\), \(P_{1}\), \(P_{2}\), and \(\omega\). Therefore, the PINN results analysis will primarily concentrate on these four states, denoted by the graphs’ blue lines. Moreover, for all cases, the training dataset consists of \(48\) data points for data loss with \(100\) collocation points for physics loss. The plots show the mean value of results for 30 realization of PINNs.
ence. We believe that the observed discrepancies between the experimental data and model predictions could be attributed to multiple factors. These include model uncertainty arising from approximations, simplifications, assumptions in the model, and the influence of noise in the experimental data.

### Parameter estimation results

The results for parameter estimation in the experimental investigations \(1\) and \(2\) across three cases are presented in the following sections. It should be noted that we did not measure the emulsion bulk modulus during the experiments. As a result, we cannot assess the accuracy of the values determined by the PINN in the experimental data scenario. The viscosity values obtained were compared to the Brinkman [44] model, described in Section 2.2. For density estimation, we relied on measurements from the Coriolis meter. In this section, we have discussed absolute percentage error, but additional tables presenting actual values can be found in Section 4.3.

In this section, the mean of the absolute percentage error (MAPE) for the parameter estimation results is defined as

\[\text{MAPE}_{\Lambda}=\frac{1}{N^{init}}\sum_{i=1}^{N^{init}}\left|\frac{ \Lambda-\hat{\Lambda}_{i}}{\Lambda}\right|, \tag{20}\]

where \(N^{init}\) is the number of PINN random initializations, considered \(30\) in this study. The variable \(\Lambda\) represents the true value of the unknown parameter, and \(\hat{\Lambda}_{i}\) is the estimated unknown parameter for the \(i\)-th PINN initialization. For simplicity, we refer to the error in this section as the absolute percentage error.

#### 4.2.1 Case 1: Three unknown parameters \(B\), \(\mu\) and \(\rho\) for different types of data

In this section, we focus on Case 1, where we considered the bulk modulus (\(B\)), viscosity (\(\mu\)), and density (\(\rho\)) as unknown while considering the remaining parameters as known. In Table 8, we present the MAPE for the two experimental investigations and all data scenarios.

As can be seen in Table 8, the PINN could successfully estimate the flow parameters for all cases. In the experimental data scenario, the MAPE were slightly higher, for investigation \(1\) the errors were below \(2\,\%\), while investigation \(2\) had errors below \(4.5\,\%\). It is important to mention that for the experimental cases, the error refers to the difference between the obtained viscosity and the Brinkman [44] model, which may be not entirely accurate. Also, we considered the single effective viscosity assumption, as described in Section 2.2. Thus, this higher percentage difference on the can also be related to these assumptions and limitations.

In addition to the mean analysis, it is important to assess the parameter absolute percentage error dispersion. The error dispersion for investigation \(1\) and \(2\) is visualized in Fig. 7(a) and Fig. 7(b), respectively.

In Fig. 7(a) and Fig. 7(b), the fluid density (\(\rho\)) and viscosity (\(\mu\)) parameters demonstrate consistent behavior across multiple initializations and data sources with relatively low dispersion. On the other hand, Investigation 2 presented outliers in all parameters for the simulated case, and for \(\rho\) and \(\mu\) in the simulated with noise scenario.

When comparing Investigations 1 and 2 in Table 8 simulated and simulated with noise, we can observe that the performance in Investigation 1 is superior to that in Investigation 2. However, the box plots for both investigations suggest a similar data spread. This discrepancy can be attributed to the influence of outliers on the mean value calculation. As shown in Fig. 7(b), the mean value, denoted by a triangle symbol, lies outside the interquartile range, which suggests the impact of outliers. Furthermore, when we compared the median of the errors they presented similar values, which further indicated the influence of these outliers on the mean.

\begin{table}
\begin{tabular}{l l r r r} \hline \hline  & Data type & \(B\) & \(\mu\) & \(\rho\) \\ \hline Investigation 1 & Sim. & \(0.31\,\%\) & \(0.03\,\%\) & \(0.02\,\%\) \\  & Noisy & \(0.58\,\%\) & \(0.12\,\%\) & \(0.09\,\%\) \\  & Exp. & N/A & \(1.80\,\%\) & \(0.80\,\%\) \\ \hline Investigation 2 & Sim. & \(0.61\,\%\) & \(0.32\,\%\) & \(0.29\,\%\) \\  & Noisy & \(0.48\,\%\) & \(0.21\,\%\) & \(0.21\,\%\) \\  & Exp. & N/A & \(4.30\,\%\) & \(2.40\,\%\) \\ \hline \hline \end{tabular}
\end{table}
Table 8: **Mean of the absolute percentage error for Case 1.**
Notably, the addition of noise negatively affected the results, as for most of the parameters, their distributions presented a slightly wider dispersion than the simulated. Further, analyzing the experimental data results, it is evident that they exhibit a wider distribution when compared to the simulated and simulated-with-noise data. This is expected as experimental data is more challenging due to noise and potential model inaccuracies.

2.2 Case 2: Seven unknown parameters \(B\), \(\mu\), \(\rho\), \(k_{u}\), \(k_{d}\), \(k_{3p}\) and \(k_{4p}\) for different types of data

In this section, we examine Case 2. We treated the fluid's bulk modulus \(B\), viscosity \(\mu\), and density \(\rho\), along with the pump's viscous flow loss coefficient \(k_{3p}\), the pump equivalent resistance \(k_{4p}\), and the equivalent resistances of the upstream and downstream pipeline \(k_{u}\) and \(k_{d}\) respectively, as unknowns. The remaining parameters were considered known. The MAPE from the two experimental investigations across all data scenarios is provided in Table 9.

The results presented in Table 9 demonstrate the PINN's successful estimation of flow parameters across most cases, except for the first investigation's simulated-with-noise and experimental data scenarios. In these cases, the mean of errors was significantly higher compared to the simulated, simulated and Investigation 2 simulated-with-noise and experimental data scenario. Notably, for the experimental data scenario, the viscosity (\(\mu\)) exhibited a difference of \(120.19\,\%\), and most other parameters had errors higher than \(27\,\%\). However, for the second investigation, the mean of errors in Table 9 were all below \(15.5\,\%\), indicating better performance in the experimental data scenario. As discussed earlier, the higher

Figure 8: **Box plot of the absolute percentage errors for Case 1. Each subfigure shows a specific experimental investigation and contains subplots for each unknown parameter under analysis. Within these subplots, three distinct box plots are presented: blue for the experimental data, red for the simulated data, and gray for the simulated data with added noise. These box plots include whiskers, outliers, and the mean value, represented by a triangle symbol. Outliers are identified as data points that lie more than 1.5 times the interquartile range (IQR) above the third quartile or below the first quartile. The x-axis represents the absolute percentage error, calculated with respect to the reference value. This graphical structure is consistently employed in subsequent figures related to inverse problems.**
error in the experimental case may be caused due to model uncertainty, simplifications in the model, and experimental noise.

However, the results for Investigation 1 in the simulated data scenario require some attention. As shown in Fig. 8(a), all unknown parameters presented outliers. Furthermore, their computed means were observed to lie outside the IQR. It was also previously noted in Case 1 (Section 4.2.1), which implies an influence of these outliers on the mean value. When examining the median values illustrated in Fig. 9, the error metrics for simulated scenarios in both investigations were predominantly found to be below \(1\,\%\).

The addition of noise considerably affected the results in both experimental investigations, leading to higher mean errors as shown in Table 9, whereas in investigation, the mean of the errors was smaller than in experimental Investigation 1. Analyzing the box plots in Fig. 9, we can observe that the simulated noisy scenario resulted in greater dispersion for most parameters. This observation indicates that the introduction of noise adversely affected the accuracy of the estimated parameters. Moreover, we can observe that similar to the discussed in Section 4.2.1, the experimental data scenario presented the highest mean errors across the initializations for both experimental investigations, except for the \(k_{4p}\) in the second experimental investigation.

Figure 9: **Box plot of the absolute percentage errors for Case 2. The box plot illustrates the distribution of absolute percentage errors for each parameter in Case 2. Each subplot represents an unknown parameter for the data scenarios evaluated. The x-axis has been changed to logarithmic for better visibility. This graphical structure is consistently employed in subsequent figures related to inverse problems.**
Despite the high errors in Investigation 1's experimental data scenario, most estimated parameters exhibited relatively small dispersion across the weight initializations, as evident in Fig. 9a. It is worth noting that the parameter selection for this case did not consider the practical identifiability analysis, which assessed, in this study, the parameter sensitivity to noise. Therefore, this higher mean error with relatively low dispersion suggests that the system might not be practically identifiable within the flow condition of experimental investigation 1. It is possibly converging to a parameter set other than the true one. It is important to note that as the water cuts increase up to the emulsion phase inversion, the fluid behavior becomes more complex, posing a more challenging scenario than lower water cuts. Further, in our practical identifiability analysis, we considered only the sensitivity analysis rather than the Monte Carlo approach, where the optimization algorithm, in this case, the PINN, would be also included.

Therefore, to draw more definitive conclusions about how the parameter estimation is sensitive to the water cut, additional investigations need to be conducted with a larger experimental test matrix that covers a wider range of water cuts. Additionally, a Monte Carlo approach with the PINN on the practical identifiability analysis could potentially better evaluate the effects of the model uncertainties with the experimental data scenario, as the approach described in Section 3.2 considered only the simulated data. Further, as the local structurally identifiability analysis assesses the identifiability of the neighborhood of the parameter set, it is possible that restricting the search range and imposing stricter bounds on the parameter search could lead to better results. This approach is suggested in [49].

2.3 Case 3: Eight unknown parameters \(B\), \(\mu\), \(\rho\), \(k_{u}\), \(k_{d}\), \(k_{4p}\), \(k_{1s}\) and \(k_{5s}\) for different types of data

In this section, we present the results of the inverse problem for Case 3, specifically evaluating the performance of the parameter estimation. The evaluation is based on the MAPE, as shown in Table 10. The considered parameters, in this case, include the bulk modulus (\(B\)), viscosity (\(\mu\)), and density (\(\rho\)) of the fluid, as well as the pump equivalent resistance (\(k_{4p}\)), the first impeller-fluid coupling coefficient (\(k_{1s}\)), the shaft second-order friction coefficient (\(k_{5s}\)), and the equivalent resistance of the pipeline upstream and downstream (\(k_{u}\) and \(k_{d}\)).

The results in Table 10 and Fig. 10 indicate that the PINN successfully estimated the density and viscosity for all scenarios in Case 3, but it encountered challenges in estimating the bulk modulus. The \(15\,\%\) value represents the upper and lower bounds imposed on the output scaling transformation for the bulk modulus. Additionally, when analyzing the dispersion in Fig. 10, it is evident that the distributions of the simulated with added noise case are similar to those of the simulated case, with similar IQR.

In Table 10, the distinction between data sources with smaller errors, as observed in Cases 1 and 2 (Section 4.2.1 and Section 4.2.2), is not evident. The Investigation 1 and 2 demonstrated errors mostly below \(13\,\%\). Notably, in Investigation 1, the upstream pipeline equivalent resistance (\(k_{u}\)) exhibited the highest errors across all data sources, reaching up to \(34.51\,\%\). However, in Fig. 3b, we can observe that the \(k_{u}\) presented a relatively high correlation with other parameters, such as the \(k_{d}\). On the other hand, in Investigation 2, we do not observe errors as high as in Investigation 1. This discrepancy can likely be attributed to the distinct flow characteristics of each investigation. As previously discussed in Section 4.2.2, further tests under diverse flow conditions are essential to understand the performance difference between Investigation 1 and 2 comprehensively.

It is noteworthy that, for Case 3, besides the fluid bulk modulus and density, the other parameters also have bounded output scaling transformations, which could have contributed to lower errors. Thus, although Case 3 poses an additional challenge compared to Case 2, once it involves estimating eight unknown parameters instead of seven, Case 3 and Case 2 are not directly comparable.

\begin{table}
\begin{tabular}{l l r r r r r r r r} \hline  & & \(B\) & \(\mu\) & \(\rho\) & \(k_{u}\) & \(k_{d}\) & \(k_{4p}\) & \(k_{1s}\) & \(k_{5s}\) \\ \hline Inv. 1 & Sim. & \(14.95\,\%\) & \(11.31\,\%\) & \(3.02\,\%\) & \(26.27\,\%\) & \(3.25\,\%\) & \(1.79\,\%\) & \(6.92\,\%\) & \(12.24\,\%\) \\  & Noisy & \(14.51\,\%\) & \(10.41\,\%\) & \(1.88\,\%\) & \(25.97\,\%\) & \(4.99\,\%\) & \(2.06\,\%\) & \(6.88\,\%\) & \(6.00\,\%\) \\  & Exp. & N/A & \(13.05\,\%\) & \(9.26\,\%\) & \(34.51\,\%\) & \(16.01\,\%\) & \(5.87\,\%\) & \(12.37\,\%\) & \(7.76\,\%\) \\ \hline Inv. 2 & Sim. & \(15.00\,\%\) & \(14.34\,\%\) & \(4.12\,\%\) & \(6.90\,\%\) & \(6.66\,\%\) & \(3.43\,\%\) & \(12.71\,\%\) & \(6.96\,\%\) \\  & Noisy & \(15.00\,\%\) & \(13.66\,\%\) & \(3.77\,\%\) & \(4.92\,\%\) & \(13.63\,\%\) & \(7.04\,\%\) & \(11.86\,\%\) & \(12.22\,\%\) \\  & Exp. & N/A & \(1.60\,\%\) & \(10.16\,\%\) & \(1.96\,\%\) & \(12.05\,\%\) & \(5.28\,\%\) & \(7.01\,\%\) & \(9.25\,\%\) \\ \hline \end{tabular}
\end{table}
Table 10: **Mean absolute percentage error for Case 3.**
### Estimated unknown parameters mean and standard deviation

In this section, we present additional PINN results for parameter estimation in the ESP system model. These results are obtained by averaging parameter estimations over the PINN 30 training runs and computing the standard deviation. The tables Tables 11 and 12 contain the results of the first and second experimental investigations, respectively. It is important to note that these tables serve as complementary results of the discussed in Section 4.2 and are not the focus of detailed discussion in this work. They are included for the sake of completeness and may be of interest to future research works. Any parameter not accounted for in a specific case is denoted by the symbol.

For each of the three cases, defined based on varying assumptions regarding known parameters, the table presents the performance of PINN in estimating the unknown parameters across the three data scenarios. We average each parameter estimation over the 30 weight initializations, and the standard deviation is indicated by the symbol \(\pm\). Parameters not considered in a given case are denoted by the \(-\) symbol.

Figure 10: **Box plot of the absolute percentage errors for Case 3. The box plot illustrates the distribution of absolute percentage errors of each parameter in Case 3. The individual subplots correspond to the unknown parameters across the evaluated data scenarios. To enhance visibility, the x-axis has been changed to a logarithmic scale.**
## References
Summary

In this study, our primary focus was on developing a physics-informed neural network (PINN) model in order to estimate unknown states and parameters of the electrical submersible pump (ESP) system operating in a two-phase flow condition when only limited state variables are known. In the present study, we only know the pump intake and discharge pressure states measurements. The PINN model may reduce the field laboratory test required to estimate the fluid properties. We first investigated the structural and practical identifiability of the ESP system, defining three distinct sets of unknown parameters. Subsequently, we assessed the performance of the PINN model on these parameter sets across two experimental studies (water to oil ratio) encompassing three data scenarios: simulated, simulated with added noise, and experimental data.

We conducted a local structural identifiability analysis using the differential elimination method for dynamical models via projections, as proposed by Dong _et al._[23]. The results revealed that the system is not locally structurally identifiable when both intake and discharge pressures are the only known parameters. For local structural identifiability, we had to consider the parameters like shaft inertia (\(I_{s}\)), fluid-impeller disk friction (\(k_{3s}\)), and shaft viscous damping (\(k_{4s}\)) as known, resulting in \(12\) identifiable unknown parameters. Then, we derived the parameter correlation matrix using the Fisher information matrix for the practical identifiability analysis. Practical identifiability provides a qualitative analysis for the unknown parameters in the presence of noise. To attain this identifiability, it became necessary to define more parameters as known. Thus, we selected the most critical parameters in the context of the ESP's actual operation as unknown and other parameters as known. These unknown parameters included the fluid properties and pipeline equivalent resistance.

The proposed PINN model in this study successfully estimated the ESP's fluid properties (bulk modulus, density, and viscosity) and dynamics of the states with reasonable accuracy. In the case when we considered fluid properties, pipeline equivalent resistances, shaft, and pump parameters, the model provides a reasonable accuracy when considering simulated data with low water cut. However, we observed a relatively higher error for higher water cut or when we considered simulated with noise and experimental data with lower water cut.

In this case, we observed outliers in the predicted results. Further study may include a robust PINN method to eliminate the outliers and better performance. Furthermore, we also observed that errors are higher in the case of higher water cuts. This error is higher in the case of noise and experimental data than in the simulated cases. These cases with more unknowns may require additional measurements of the state variables, which are generally not measured in the field. Future studies may also include further investigations with a broader range of water fractions to understand better the water cut influences. The study may also focus on updating the numerical model for higher water cuts if necessary, as assumptions made on the numerical model may not be valid for higher water cuts.

Obtaining the estimated parameters and states for the ESP system in oil production can be challenging, and the PINN offers a promising and cost-efficient alternative for estimating them. Although this approach has its benefits, it also has some limitations. One of the limitations is that the technique is computationally intensive. This means that for every desired estimation of properties, the PINN must undergo training. More efficient PINN algorithms need to be investigated for different operational conditions. Furthermore, the accuracy of the estimated properties using PINN is heavily reliant on the accuracy of the measured states and known parameters. In the event of a faulty reading or instrument failure, the PINN's ability to provide accurate estimations is compromised.

Additionally, a Monte Carlo approach with the PINN on the practical identifiability analysis could potentially better evaluate the effects of the model uncertainties with the experimental data scenario. As the experimental data contains noise, further study may also include uncertainty quantification study. Furthermore, an ESP system model that considers the multistage separately would be relevant for ESPs in actual production, whose number of stages are considerably higher than the one considered in this study. Moreover, the current model does not consider the electrical domain, and its inclusion would be more realistic. Another notable aspect is the single viscosity assumption for the system, highlighted in Section 2.2, which remains a limitation.

### Acknowledgment

We gratefully acknowledge the support from the Energy Production Innovation Center (EPIC) at the University of Campinas (UNICAMP), financially backed by the Sao Paulo Research Foundation (FAPESP, grants 2017/15736-3 and 2019/14597-5). We also acknowledge the financial sponsorship from Equinor Brazil and the regulatory support offered by Brazil's National Oil, Natural Gas, and Bio
fuels Agency (ANP) through the R&D levy regulation. Additional acknowledgments are extended to the Center for Energy and Petroleum Studies (CEPTRO), the School of Mechanical Engineering (FEM), and the Artificial Lift and Flow Assurance Research Group (ALFA).

The work of K. Nath and G. E. Karniadakis was supported by OSD/AFOSR MURI grant FA9550-20-1-0358.
