Two-Stage Wiener-Physically-Informed-Neural-Network (W-PINN) AI Methodology for Highly Dynamic and Highly Complex Static Processes

1. Introduction

Our objective is the extension of the Wiener-Physically-Informed-Neural-Network (W-PINN) approach developed by [1] to improve modeling effectiveness when processes (i.e., systems) are highly nonlinear. W-PINN is a two-step modeling methodology. In Step 1, p inputs enter and exit their own linear dynamic modeling structure. In Step 2, the p dynamic outputs enter a static modeling structure that outputs the fitted response. W-PINN effectively and accurately modeled a real industrial grain dryer, a batch process with first order (FO) dynamic behavior and simple (i.e., one input) linear regression static behavior [1]. [2] recently introduced the Theoretically-based Dynamic Regression (TDR) modeling approach to tie it to, and to distinguish it from, Dynamic Regression (DR) [3,4,5,6,7], a familiar lagged-based approach in several academic disciplines. More specifically, the core difference between TDR and DR is that TDR uses physically-based dynamic structures and DR uses lagged (i.e., empirically)-based structures. Thus, TDR is more parsimonious, physically interpretable, and exploitable than DR. However, with its limitation to first- and second-order linear dynamic structures, and its lack of connection to ANN, TDR can be viewed as a subclass of W-PINN.

[2] applied TDR in three (3) modeling cases. The first one modeled a subject’s weight (y) with the consumption of four nutrients (x_i), i.e., four inputs. The fitted model is a first-order (FO), linear dynamic/quadratic static structure. The training r_fit = 0.92 and the validation r_fit = 0.91. The second case in [2] modeled the top tray temperature (y) of a pilot distillation column with nine (9) inputs (x_i). It resulted in a complex second-order linear dynamic structure for each input and a simple FO multiple-input static regression function. The training r_fit = 0.96, the validation r_fit = 0.97, and the average testing r_fit for nine (9) data sets over a three-year period equaled 0.84. The third case modeled the sensor glucose concentration (SGC) (y) of eleven (11) type 1 diabetes data sets using twelve (12) inputs in outpatient (free-living) data collection of two weeks. It also resulted in a complex second-order linear dynamic model structure for each input. However, the static structure is quite complex and limited this case to a first-order, multiple-input, static structure. For the second study, there was one process, the distillation column. In the SGC modeling study, each subject is a process, or system, and thus has its own one-week of training/one-week of validation results. For these eleven sets of results, the average $r_{fit},$ ${\overline{r}}_{fit},$ for training and validation equaled 0.62 and 0.68, respectively.

Moreover, the objective of this work is the development of an effective W-PINN modeling approach for dynamic systems with highly complex static behavior. More specifically, using the Stage 1 SGC dynamic modeling results obtained by Rollins et al. (2025) [2], the objective of this work is to significantly increase ${\overline{r}}_{fit},$ using a novel, proposed, two-stage, W-PINN modeling approach.

The classical ANN approach is a one-box, empirical modeling methodology as illustrated in Fig. 1a. As shown, p measured x_i inputs enter the ANN function, $\hat{f} (\mathit{X})$, an estimated (as denoted by “^”) empirical function with constant coefficients that are adjusted under some criterion, commonly least squares estimation, to maximize agreement (i.e., fit) between its modeled output, $\hat{y}\hspace{0.33em}=\hspace{0.33em}\hat{f}(\mathit{X})$ and its measured output, y (e.g., SGC), where $\mathit{X}$ is a p^th dimensional vector of measured inputs. In Fig. 1a, $\hat{f} (\mathit{X})$ can be a static function, e.g., a non-linear regression function, or a combined static and empirically dynamic function of lag variables, e.g., a Long- and Short-Term Memory (LSTM) [8, 9] function.

Phenomenologically dynamic and empirically static methodologies are illustrated in Figs. 1b (PINN [10]) and 1c (W-PINN, [1]). [10] named their methodology “physically-informed-neural-network” (PINN). The critical difference between Fig. 1a (classical ANN) and Figs. 1b and 1c are the number of stages. Figure 1a has one stage for static and dynamic structures, and the two-stage methods in Figs. 1b and 1c have one stage for static structures and one stage for dynamic structures. The PINN dynamic block is not restricted to linear dynamic structures as it is for W-PINN. However, nonlinear dynamic behavior is modeled when the dynamic outputs of W-PINN (i.e., ${\hat{v}}_i\prime s$) are passed through $\hat{f} (\widehat{\mathit{V}})$. A critical advantage of W-PINN over PINN is that each input has its own dynamic model structure, as illustrated by comparing Figs. 1b and 1c.

The next section, Section 2, describes the W-PINN methodology fundamentally and mathematically. The theoretical structure of a general and complete second-order dynamic structure is first given in Section 2.1 as differential equation and then transformed to its discrete-time version using backward difference derivatives. It then gives the explicit equation for ${\hat{v}}_{i,t}^{\prime }$, where “t” is the sampling time. Section 2.2 gives important Stage 1 details and results to assist in the understanding of the Stage 2 methodology. Section 2.3 describes the JMP^® and Python^® Stage 2 methodologies. Section 2.4 give the mathematical details of Model 1 (input only model), Model 2 (input-output) and Model 1-2, a combination of the strengths of Models 1 and 2. Section 2.5 gives the forecast model structure, i.e., the structure for ${\widehat{\mathit{y}}}_{\mathit{t}+\mathit{k}\Delta \mathit{t}}.$ Finally, Section 2.6 gives the information and equations for the summary statistics.

Section 3 gives a table with all the numerical results for the three modeling methods. It also gives Model 1, Model 2 and Model 1-2 graphical results for the best fitting Stage 2 subject, Subject 2. For all three models, their validation $r_{fit}$ is 0.93. Section 4 gives a discussion of the results and Section 5 comments on work in progress and speculates on other possible future directions.

2. Materials and Methods

This section describes the two-step W-PINN methodology in detail. It also gives critical Stage 1 SGC modeling particulars used by [2] to obtain the ${\hat{v}}_{\mathit{i}\mathit{j}}\prime s$ posted on the website of the last author at https://drollins9.wixsite.com/derrickrollins. With both the $x_{ij}\prime s$ and ${\hat{v}}_{\mathit{i}\mathit{j}}\prime s$ posted on this website, modelers have the option of using these data sets to build both stages or just the second stage, the aim of this work.

With the sampling rate, Δt, equal to 5 minutes, our SGC models are forecasting 12-steps (i.e., 60 minutes) into the future. The forecast nature of this work is an artifact of the data sets we are using and their application and, thus, not a necessity of the methodology. As described in [2], 12Δt is the estimated observable time it took for the manipulated variable (MV), exogenous insulin, to cause SGC to start decreasing after a bolus increase (or insulin injection). The 60-minute estimate was very consistent for the subjects in this clinical study as noted in [2].

The models that this work develops are for an unobservant, 60-minute, forecast monitoring scenario. “Unobservant” is meant to convey the protocol that the person(s) determining insulin changes have no knowledge of the forecast (${\widehat{y}}_{ij}$) estimates. In addition, even though the authors of [2] worked diligently to obtain models that minimize ${\widehat{y}}_{ij}$ pairwise correlation, we note that it is still significantly present. Thus, for this reason, this work is best understood as an unobservant monitoring application and not applicable to closed-loop control.

2.1. W-PINN

Our W-PINN approach uses backward difference derivatives (BDD) to discretize second-order-plus-dead-time-plus-lead (SOPDTPL) (see [12]) theoretical dynamic systems, the only type used in this work, as given in Eq. 1, below. The dynamic system does not have to be initially at a steady state for our W-PINN modeling methodology since the initial conditions are also estimated. Note that, Eq. 1 is the expression for each of the p-inputs.

$${\tau }_i^2{\displaystyle \frac{d^2{v}_i(t)}{d{t}^2}}+2{\tau }_i{\zeta }_i{\displaystyle \frac{d{v}_i(t)}{dt}}+v_i(t)={\tau }_{ai}{\displaystyle \frac{d{x}_i(t-{\theta }_i)}{dt}}+x_i(t-{\theta }_i)$$

(1)

with

$$E[y(t)]=f(\mathit{V}(t))$$

(2)

where $t\ge {\theta }_i\ge 0,$ ${\tau }_i>0,$ and ${\zeta }_i\ge 0$for i = 1, …, p, x_i(t) is the value of the i^th input variable at t, and v_i(t) is the value of the i^th output variable at t, in the units of x_i, y(t) is the output variable in its units at t, $E[y(t)]$ means the expected value (i.e., true mean) of y(t), $f(\mathit{V}(t))$ is the true output (gain) function of V(t), the vector of the v_i(t)’s. When $f(\mathit{V}(t))$ is a nonlinear function of V(t), as in the ANN (i.e., W-PINN) case, Eqs. 1 and 2, taken together, have a Wiener block-oriented structure [11], as shown in Fig. 1c.

The lead term is the first term on the right side of the equal sign in Eq.1. This term tends to "speed up" the response and provides what the process modeling and control community has termed "numerator dynamics" [11,12,13]. [11] developed a second-order, multiple-input, single-output, discrete-time, nonlinear Wiener dynamic approach using BDD based on Eq. 1. More specifically, using BDD approximation applied to a sampling interval of Δt, an approximate discrete-time form of Eq. 1 is

$${\widehat{v}}_{i,t}=\left\{\begin{array}{l}{\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t-\Delta t}+{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t-2\Delta t}+\hspace{0.33em}{\widehat{\omega }}_{1,i}{x}_{i,t-{\widehat{\theta }}_i-\Delta t}+{\widehat{\omega }}_{2,i}{x}_{i,t-{\widehat{\theta }}_i-2\Delta t},\hspace{1em}\hspace{1em}t>{\widehat{\theta }}_i=\widehat{m}{}_i\Delta t\\ {\widehat{v}}_{{\widehat{\theta }}_i}={\widehat{v}}_{{\widehat{m}}_i\Delta t}\hspace{0.33em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{1em}\hspace{0.33em}t={\widehat{\theta }}_i={\widehat{m}}_i\Delta t\\ \hspace{1em}\hspace{0.33em}is\hspace{0.33em}undefined,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.33em}\hspace{0.33em}t<{\widehat{\theta }}_i={\widehat{m}}_i\Delta t\end{array}\right.$$

(3)

with

$${\widehat{\delta }}_{1,i}={\displaystyle \frac{2{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(4)

$${\widehat{\delta }}_{2,i}={\displaystyle \frac{-{\widehat{\tau }}_i^2}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(5)

$${\widehat{\omega }}_{1,i}={\displaystyle \frac{({\widehat{\tau }}_{ai}+\Delta t)\Delta t}{{\widehat{\tau }}_i^2+2{\widehat{\tau }}_i{\widehat{\zeta }}_i\Delta t+\Delta t^2}}$$

(6)

such that ${\widehat{\omega }}_{2,i}=1-{\widehat{\delta }}_{1,i}-{\widehat{\delta }}_{2,i}-{\widehat{\omega }}_{1,i}$to satisfy the unity gain constraint. From Eq. 3 with $t>{\widehat{\theta }}_i+2\Delta t,$

$$\begin{array}{l}{\widehat{v}}_{i,t}-{\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t-\Delta t}-{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t-2\Delta t}\\ \hspace{1em}\hspace{1em}={\widehat{\omega }}_{1,i}{x}_{i,t-\left({\widehat{\theta }}_i+\Delta t\right)}+{\widehat{\omega }}_{2,i}{x}_{i,t-\left({\widehat{\theta }}_i+2\Delta t\right)}\end{array}$$

(7)

$$\begin{array}{l}\Rightarrow \hspace{1em}\left(1-{\widehat{\delta }}_{1,i}B-{\widehat{\delta }}_{2,i}{B}^2\right){\widehat{v}}_{i,t}\\ \hspace{1em}\hspace{1em}=\left({\widehat{\omega }}_{1,i}{B}^{{\widehat{\theta }}_i+1}+{\widehat{\omega }}_{2,i}{B}^{{\widehat{\theta }}_i+2}\right)x_{i,t}\end{array}$$

(8)

$$\Rightarrow \hspace{1em}{G}_{i,t}={\displaystyle \frac{{\widehat{v}}_{i,t}}{x_{i,t}}}={\displaystyle \frac{{\widehat{\omega }}_{1,i}{B}^{{\widehat{\theta }}_i+1}+{\widehat{\omega }}_{2,i}{B}^{{\widehat{\theta }}_i+2}}{1-{\widehat{\delta }}_{1,i}B-{\widehat{\delta }}_{2,i}{B}^2}}$$

(9)

After obtaining ${\widehat{v}}_{i,t}$ for each input i, the modeled output value, at time t, is determined by entering these results into $\widehat{f}({\hat{\mathit{V}}}_t)$, a static ANN in this application, i.e.,

$${\widehat{y}}_t=\widehat{f}\left({\widehat{\mathit{V}}}_t\right)$$

(10)

2.2. Stage 1 Modeling Method

This subsection gives important Stage 1 details and results to assist in the understanding of the Stage 2 methodology. While missing output (i.e., SGC) measurements are acceptable, missing input values are not for discrete-time modeling. Activity tracker data were the only missing input data. These missing values were estimated by averaging the two values on both sides of a gap and filling in the gap with this value. Some gaps were several hours long. Cross-validation [6] was used to guard against overfitting, with the first week as the training (Tr) data set, and the second week as the validation (Val) data set.

In Stage 1, all inputs were first modeled separately on their own Excel^® worksheet with a first-order linear regression static function. For each case, insulin was modeled first. The estimated deadtime, i.e., ${\widehat{\theta }}_{MV},$ was set at 60 minutes and was varied one Δt forwards and backwards at a time to find the value that gave the best fit. For all the inputs, the estimate of θ_MV, ${\widehat{\theta }}_{MV},$ was determined to be 60 minutes, i.e., 12 Δt. The food variables were the only ones with announcements, and the carbohydrate input was the only one found to have a deadtime less than ${\widehat{\theta }}_{MV}$. The “time of day” input has no deadtime, and the dead time for all other inputs was ${\widehat{\theta }}_{MV},$ except for fats that had deadtimes that were much larger than ${\widehat{\theta }}_{MV},$ as determined by model estimation. After determining the dynamic coefficient estimates for each input (i.e., Eqs. 4 to 6), these values were copied to an Excel^® worksheet as the dynamic structure starting values for fitting the SOPDTPL dynamic, and first-order static, multiple-input model (i.e., Eq. 10). With $\widehat{f}\left(\widehat{\mathit{V}}\right)$ as a first-order multiple linear regression static function, $\widehat{\mathit{V}},$ for Stage 2 was determined.

2.3. Stage 2 Model Development Modeling Methods

The objective of this work is the development, evaluation, and comparison of two, Stage 2, W-PINN modeling approaches using ${\hat{\mathit{V}}}_t$ from Stage 1 to obtain $f({\hat{\mathit{V}}}_t)$ for each of the eleven data sets for the two approaches. The first approach uses the JMP^® ANN toolbox to approximately find the smallest SSE (i.e., SSR) by fitting many cases and selecting the best one. The second approach is our newly developed, confidential, ANN structure that is coded using Python^®. As Section 3 will show, both approaches significantly improved fit in comparison to $f({\hat{\mathit{V}}}_t)$ as a first-order linear regression function.

2.4. Three Input Models

We developed three types of input model structures for this application. The first one we call the “input only model” or “Model 1.” All the inputs in this structure have a deadtime $\ge \hspace{0.33em}{\widehat{\theta }}_{MV}$ except for announcement inputs that can have deadtimes less than$\hspace{0.33em}{\widehat{\theta }}_{MV}$ like carbohydrates, equal to $\hspace{0.33em}{\widehat{\theta }}_{MV}$ like proteins, greater than $\hspace{0.33em}{\widehat{\theta }}_{MV}$ like fats, and zero like time of day, as mentioned above.

The second one we call the “input-output model” or “Model 2.” It combines the input-only structure of Model 1 (i.e., Eq. 10) with a model of weighted residuals, a minimum of $\hspace{0.33em}{\widehat{\theta }}_{MV}$ distance in the past (note that, this is model building and not model forecasting), as shown in Eq. 11 below (see [11] for the derivation).

$$\begin{array}{l}{\widehat{y}}_t=f\left({\hat{\mathit{V}}}_t\right)+{\widehat{\varphi }}_1\left(y_{t-{\widehat{\theta }}_{MV}}-{\widehat{y}}_{t-{\widehat{\theta }}_{MV}}\right)\\ \hspace{1em}\hspace{1em}+{\widehat{\varphi }}_2\left(y_{t-{\widehat{\theta }}_{MV}-\Delta t}-{\widehat{y}}_{t-{\widehat{\theta }}_{MV}-\Delta t}\right)+\dots \\ \hspace{1em}\hspace{1em}=f\left({\widehat{\mathit{V}}}_t\right)+{\widehat{\varphi }}_1{e}_{t-{\widehat{\theta }}_{MV}}+{\widehat{\varphi }}_2{e}_{t-{\widehat{\theta }}_{MV}-\Delta t}+\dots \end{array}$$

(11)

Equation 11 has no value if any residual is not determinable due to missing output measurements. Thus, unlike Model 1, which has estimates for all t since it uses only input data, Model 2 will not have an estimate when an output value is missing.

The final model, Model 1-2, is a combination of the strengths of Models 1 and 2. More specifically, for Model 1-2,

$${\widehat{y}}_t=\left\{\begin{array}{l}\mathrm{Model} 2 {\widehat{y}}_t, \mathrm{if} \mathrm{its} \mathrm{value} \mathrm{exists}\\ \mathrm{Model} 1 {\widehat{y}}_t, \mathrm{if} \mathrm{its} \mathrm{Model} 2 \mathrm{value} \mathrm{does} \mathrm{not} \mathrm{exists}\end{array}\right.$$

(12)

2.5. Forecast Structures

Equation 10, ${\widehat{y}}_t=f\left({\widehat{\mathit{V}}}_t\right)$, is the fitted structure for Model 1, i.e., estimate of the output, $y_t$, at the current time, t. There are no missing input values in Eq. 10. This is why missing armband data had to be estimated. In addition, non-announcement input values to obtain Eq. 10 must be at least a distance of ${\widehat{\theta }}_{MV}$ in the past. This requirement is because the model developed input lag must be the same as the forecast input lag, i.e., $\widehat{\mathit{V}}$ at t, for forecasting a ${\widehat{\theta }}_{MV}$ distance into the future.

After obtaining ${\widehat{\mathit{V}}}_t$, its transformation into the kΔt forecast form, i.e., the online version, is given by Eq. 13 below:

$${\widehat{y}}_{t+k\Delta t}=f\left({\widehat{\mathit{V}}}_{t+k\Delta t}\right)$$

(13)

where, from Eq. 7, with t = t + kΔt,

$$\begin{array}{l}{\widehat{v}}_{i,t+k\Delta t}={\widehat{\delta }}_{1,i}{\widehat{v}}_{i,t+(k-1)\Delta t}+{\widehat{\delta }}_{2,i}{\widehat{v}}_{i,t+(k-2)\Delta t}\\ \hspace{1em}\hspace{1em}\hspace{1em}+{\widehat{\omega }}_{1,i}{x}_{i,t-\left({\widehat{\theta }}_i-(k-1)\Delta t\right)}+{\widehat{\omega }}_{2,i}{x}_{i,t-\left({\widehat{\theta }}_i-(k-2)\Delta t\right)}\end{array}$$

(14)

Note that, if k = 12 and Δt = 5 minutes, the Eq. 14 model is forecasting 60 minutes into the future. Thus, all the non-announcement inputs must have, i.e., use, a model building and forecast prediction dead time of at least 60 minutes.

2.6. Statistical Analyses

The first, and most important, modeling statistic is r_fit (which is bounded between -1 and 1), the fitted correlation of the measured SGC, $y_i$, and the fitted SGC, ${\widehat{y}}_i$, and given in Eq. 15 below.

$$r_{fit}=r_{y_t,{\widehat{y}}_t}={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left(y_i-\overline{y}\right)\left({\widehat{y}}_i-\overline{\widehat{y}}\right)}}{\sqrt{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{y}\right)}^2}}\cdot \sqrt{{\displaystyle \sum _{i=1}^n{\left({\widehat{y}}_i-\overline{\widehat{y}}\right)}^2}}}}$$

(15)

where n is the number of samples in the set and the bar above a statistic means that it is its sample mean value. The equations to determine AAD and AD are, respectively,

$$AAD={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left|y_i-{\widehat{y}}_i\right|}}{n}}$$

(16)

$$AD={\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left(y_i-{\widehat{y}}_i\right)}}{n}}$$

(17)

The equation for SSE (i.e., SSR, the sum of squared residuals), the more common name and used by JMP^®, is

$$SSE=SSR={\displaystyle \sum _{i=1}^n{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(18)

3. Results

Training and Validation Stage 1 and Stage 2 summary statistics for the eleven subjects are given in Table 1. All the results are r_fit unless indicated otherwise. Recalling that, each subject has a fixed ${\widehat{\mathit{V}}}_t$ that was determined in Stage 1 using the first-order static structure as given in Eq. 19, below:

$${\widehat{y}}_t={\widehat{a}}_0+{\widehat{a}}_1{\widehat{v}}_{1,t}+\dots +{\widehat{a}}_p{\widehat{v}}_{p,t}$$

(19)

As shown in Table 1, Stage 1, Model 1, r_fit,val results varied from 0.59 to 0.77, with a mean of 0.68. Moreover, Stage 2, Model 1, r_fit,val results improved significantly over the Stage 1 results for both ANN approaches. As shown, JMP^® Stage 2, Model 1, r_fit,val results varied from 0.60 to 0.85, with a mean of 0.74. However, Python^® Stage 2, Model 1, r_fit,val results are significantly better than JMP^®, varying from 0.72 to 0.93, with a mean of 0.82. As a result, Model 2 training and validation results, and Models 1-2 validation results are given in Table 1 for Python^® only. From Model 1 to Model 2, the Python^® mean r_fit,val increased from 0.82 to 0.87, the minimum from 0.72 to 0.80, and the maximum of 0.93 did not change. In summary, Python^® Stage 2 results improved considerably over Stage 1 results and are significantly better than JMP^® Stage 2 results.

Graphical Python^® Stage 2 fitted and measured SGC results for Subject 2 (the best case) are given (i.e., plotted) in Figs. 2 to 4. Figure 2 is Model 1 training and validation. Figure 3 is Model 2 training and validation. Figure 4 plots are the combined validation results, with Model 1 plotted when there is no output data, and Model 2 is plotted when there is output data, i.e., the Model 1-2 validation plot. The Model 1-2 plots are associated with the results in the last three columns in Table 1. Figure 3 shows excellent fit of Model 2 and the highly realistic behavior of Model 1 when Model 2 results are not possible because of missing SGC data (see Eq. 11).

4. Discussion

In this study, two physically-based virtual forecasting sensor approaches were developed for obtaining the value of the response variable, SGC, a θ_MV time distance in the future, for a two-stage forecast modeling application. The first stage, a physically (i.e., theoretical) based dynamic modeling approach [11] estimates the physically interpretable dynamic parameters from the measured inputs (x_i’s) with multiple physical constraints to obtain dynamic outputs (v_i’s). The v_i’s are the inputs to the second stage, a static ANN structure. For the first method, this structure was determined by using the ANN toolbox in JMP^®. For the second method, this structure was determined by using a confidential method that this work developed and coded using Python^®. Both methods resulted in large average improvements over the Stage 1 results using a first-order linear regression static structure (see Table 1). In addition, a critical advantage of these two approaches is that the modeling is much easier and much less time-consuming than the 2^nd order multiple linear regression (MLR) approach. We do note, however, that static behavior of these data sets is highly nonlinear. Thus, we strongly recommend ANN over MLR for the static model structure, i.e., ${\widehat{y}}_t=\widehat{f}\left({\widehat{\mathit{V}}}_t\right).$ Note that ANN modeling is just a particular class of multiple nonlinear regression. We also note that the MLR model was applied to Subject 11 to compare the performance with the ANN. The r_fit,val had a modest improvement from Stage 1 alone, going from 0.74 to 0.79, but much less than the 0.85 obtained by P-ANN.

We were pleasantly surprised by the P-ANN achievements of Models 1, 2, and 1-2. Model 1 has two subjects over the r_fit,val goal of 0.90 and a mean of 0.82. Model 2 has four subjects meeting or exceeding the r_fit,val goal and a significant increase in the mean r_fit,val of 0.87. The combined Model 1 and Model 2 approach, i.e., Models 1-2, had essentially the same summary statistics results as Model 2, as shown in Table 1. Thus, combining Models 1 and 2, to have continuous forecasting without missing fits did not adversely affect r_fit,val relative to Model 2, which had missing fits due to missing SGC measurements.

5. Conclusions and Future Work

The W-PINN approach is particularly powerful because each input x_i is dynamically transformed to its v_i counterpart and is the input to a static ANN. The proposed two-stage W-PINN approach greatly improved the SGC model fit for eleven historical diabetes data sets. A one-stage W-PINN approach, in its evaluation stage, is the next step in this research.

During his time as a professor, the second author gained valuable insight into the limitations of empirical modeling through a real-world industrial application. A BS Chemical Engineering student, also pursuing an MS in Statistics, undertook a summer project at a leading Midwest chemical company, which was approved for her MS thesis. The project focused on developing a multivariate Statistical Process Control (SPC) monitoring methodology for a process line. Data were collected, and an SPC control chart was developed, resulting in an excellent model fit. However, when the process exceeded control limits, adjustments to the manipulated variable based on this model failed to restore control. A subsequent attempt with new data and a revised control chart, despite another excellent fit, similarly failed to correct deviations when applied in a feedback control scenario. This experience highlighted that the control chart, designed for monitoring, was unsuitable for feedback control due to its reliance on empirical correlation rather than cause-and-effect relationships. Empirical SGC modeling, which uses free-living data and non-physiological structures, faces similar limitations, as it cannot adequately capture cause-and-effect dynamics critical for model-based control applications like automatic forecast control. In contrast, physically-informed modeling, which integrates physiological information and structure with free-living data, offers inherent intelligence and robust structure for potentially developing effective models for control applications. Our W-PINN methodology proposed in this manuscript is approaching the goal for the diabetes data set, but more research and creative screening of inputs are needed to fully realize the goal for SGC closed-loop feedback control. In addition, this advancement relies on a dual hormone scenario; one to decrease SGC, i.e., insulin; and one to as effectively and safely increase SGC, possibly glucagon. Nonetheless, our proposed two-stage methodology has promise, it seems, for physically-based, highly nonlinear static, systems or processes.

Type 1 diabetes SGC modeling for monitoring can be effective (i.e., informative) using empirical or physically-informed dynamic modeling approaches. Closed-loop Type 1 SGC automatic control is inherently forecast automatic control because a change in MV, injected insulin, will take a time of θ_MV to start lowering SGC. For automatic closed-loop control, empirical dynamic modeling approaches are not likely to succeed because they lack a cause-and-effect relationship, unlike PINN approaches. Insulin is the process variable that is changed to keep SGC close to its set point, i.e., it is the manipulated variable (MV). For a control system to do this well in a forecast feedback control scheme, the controlled variable, $S{G}_{\theta MV},$ must be accurately estimated. An empirical method could possibly control $S{G}_{\theta MV}$ online accurately if the correlation structure remains the same as it was when the model was developed. However, it is not possible for the correlation structure of an empirical forecast modeling approach in this context to remain intact, i.e., fixed, in online forecast feedback control because the correlation structure changes each time the controller signal to the manipulated variable is transmitted. Thus, it is prudent to restrict free-living empirical modeling to monitoring open-loop processes but not to make decisions on how much to change a manipulated variable to make changes in the control variable.

There are several challenges of the data sets used in this work. First, they are nearly a decade and a half old, and SGC technology has advanced considerably, particularly in terms of missing and lost data. Secondly, wearable technology have advanced considerably in reliability, measured sensor technology such as heart rate, as well as in data management. In addition, there are advancements in ways to get accurate consumption of food nutrients. Thus, one (distant) future goal is to evaluate W-PINN using data generated by current technology.