Developer Handoff — PINN for Hydrate Equilibrium (Teq) prediction

Source: OTC/SPE paper supplied (OTC-35362-MS.pdf). See Methodology / Physical Model

sections (≈ pages 5–7) for the full derivation and references (Bishnoi intrinsic kinetics + Towler &
Mokhatab correlation for Teq).
Goal: produce a robust PINN implementation that, given pipeline/state features, predicts
hydrate equilibrium temperature Teq(P,features) (or equivalently predicts Teq for a sweep of
Pressure) while enforcing physics (Bishnoi/kinetics residual) and returning confidence estimates.

This document is a developer-ready specification derived from the supplied paper and the
notebook outputs. It describes the exact data formats, physics formulas to implement, the PINN
architecture, loss functions, diagnostics, tests, and deliverables required. Give this file to a
developer and they should be able to implement a robust PyTorch PINN pipeline that trains and
predicts hydrate equilibrium temperature (Teq) versus pressure.

1) Summary (layman)

Input data: measurements and operating features (pressure, temperature, time, gas gravity,
pipe geometry, hold-up, etc.) from training runs or synthetic training set.

Target (what to predict): hydrate equilibrium temperature Teq (°C or K) as a function of

operating conditions — typically Teq(P,Sg,...) or the PVTSim curve (pressure vs temperature).

Method: train a neural network that (a) fits measured Teq from data (data loss) and (b) is
constrained by known physics (the hydrate-formation kinetics / Bishnoi equation) via an
additional physics residual term in the loss. This is a Physics-Informed Neural Network (PINN).

Output: predicted Teq for a pressure sweep, plus an uncertainty (confidence) measure per
predicted point (either as predicted stddev from probabilistic head or ensemble-based std).

a) Primary curve file (the "truth" hydrate curve from PVTSim)

Filename suggested: hydrate_curve_clean_fixed.csv

Columns (required):

 Pressure — pressure at curve point, units = bara (numeric)

 Temperature — equilibrium temperature at that pressure, units = °C (numeric)

(If the PVTSim export uses Kelvin, convert to °C by subtracting 273.15 before saving.)

Example 8-row snippet (CSV):

b) Training dataset (if using dynamic/flow data)

Filename: fair_pinn_synth_data.csv (or training_data.csv)

Columns (recommended) — order is important and must match scaler/model:

 time_s (s or seconds)
 P_bar (pressure in bara)
 T_K (temperature in Kelvin) — or use T_C in °C but be consistent
 Sg (gas gravity — dimensionless)
 phase_frac (gas/liquid fraction as used in their features)
 pipe_d_m (pipe diameter meters)
 holdup_liq (liquid holdup fraction)
 DeltaT (if used in synthetic features)
 y_hh (hydrate holdup) optional depending on model
 y_dpr (DPR head target) optional
 (Any other features used in the original model — maintain the same names and order.)

c) Saved artifacts to load for inference (developer should produce)

 scaler.pkl or teq_scaler.pkl (sklearn StandardScaler or equivalent) — must contain

n_features_in_ or clearly documented feature order
  pinn_model_teq_aug.pth (PyTorch saved state_dict or saved dict with state_dict,
normalization stats)

 cleaned_data.pth (optional metadata mapping X_cols or template medians)

 Output: hydrate_curve_prediction_teq_regen.csv (predictions)

3) Physics / formulas the developer must implement

3.1 Bishnoi intrinsic kinetics (paper excerpt)

The paper uses the Bishnoi intrinsic kinetics for hydrate growth (the gas consumption rate
equation). The text indicates a form like:

dmgdt=−k1⋅AS⋅hh⋅f(T,Teq)\frac{dm_g}{dt} = -k_1 \cdot A_S \cdot h_h \cdot f(T, T_{eq})dtdmg

=−k1⋅AS⋅hh⋅f(T,Teq)

Where:

 mgm_gmg = mass of gas consumed (guest)

 ttt = time

 ASA_SAS = contact surface area between phases (may be approximated e.g. via Biberg
correlations)

 TTT = system temperature

 TeqT_{eq}Teq = hydrate equilibrium temperature (the unknown Teq we want)

 k1,k2k_1, k_2k1,k2 = kinetic constants (paper lists fitted constants — developer must
copy exact numeric values from the paper or re-fit)

 hhh_hhh = hydrate holdup (volume fraction)

The paper substitutes density assumptions and simplifies so the right-hand side becomes an
expression that depends on Teq, T, and other known constants. The PINN imposes that the
predicted Teq satisfies the residual of this equation (i.e., physics residual = 0) given measured T
and other parameters.

Implementation note: Use the exact equation from the paper’s Methodology section. The
derivation reduces to a residual formula R(P,Teq,features) that the PINN enforces via a MSE of
residuals.

3.2 Towler & Mokhatab correlation (Teq estimate)

  The paper uses the Towler & Mokhatab correlation to approximate Teq from P and gas
gravity Sg where needed. Developer should implement that correlation directly (copy
the formula and units from paper or the original Towler & Mokhatab reference).

 If that correlation is not implemented, the paper falls back to using the measured
PVTSim curve for Teq in the physics residual.

If developer cannot access the exact Towler& M formula right away: use the PVTSim curve (or
a lookup/interpolation) and use Towler for sanity checks only.

4) PINN architecture & loss (exact instructions)

4.1 Network

 Use PyTorch.

 Input dimension D = number of features (e.g. 8). The developer must maintain the same
X_cols order as the scaler.

 Hidden layers: 4 hidden layers, 128 units each (ReLU). (Paper used similar sizes —
developer can tune.)

 Output head: probabilistic Teq head with two outputs:

o mu (predicted mean Teq) — linear output

o logvar (predicted log variance) — linear output (to enforce positive variance: std
= sqrt(exp(logvar)))

(Alternatively: two outputs mu and logvar, or mu + sigma with softplus activation for sigma.)

class PINN_Teq(nn.Module):

def __init__(self, input_dim):

super().__init__()

self.body = nn.Sequential(

nn.Linear(input_dim, 128), nn.ReLU(),

nn.Linear(128, 128), nn.ReLU(),

 nn.Linear(128, 128), nn.ReLU(),

nn.Linear(128, 128), nn.ReLU(),

self.mu_head = nn.Linear(128, 1)

self.logvar_head = nn.Linear(128, 1)

def forward(self, x):

h = self.body(x)

mu = self.mu_head(h) # normalized units (if input normalized)

logvar = self.logvar_head(h)

return mu, logvar

4.2 Normalization

 Standardize inputs: StandardScaler() fit on training X (save mean/std).

 Normalize outputs (Teq): scale to zero mean & unit std for training stability (model
predicts normalized mu_norm and logvar — convert back at inference).

4.3 Loss function

Total loss = data loss + λ_phys * physics loss + λ_reg * weight decay (optional)

Data loss (negative log-likelihood for predicted Gaussian):

For each sample i:

Ldata,i=12exp⁡(−log⁡σi2) (yi−μi)2+12log⁡σi2\mathcal{L}_{\text{data},i} = \frac{1}{2}\exp(-\log \

sigma_i^2)\,(y_i - \mu_i)^2 + \frac{1}{2}\log \sigma_i^2Ldata,i=21exp(−logσi2)(yi−μi)2+21logσi2

Sum/mean over batch.

Physics loss (residual MSE):

Let R(P,mu,features) be the physics residual derived from Bishnoi equation (replace Teq with
model mu in proper units -> get residual). Then

Lphys=1N∑iR(Pi,μi,featuresi)2\mathcal{L}_{\text{phys}} = \frac{1}{N}\sum_i R(P_i, \mu_i, \

text{features}_i)^2Lphys=N1i∑R(Pi,μi,featuresi)2
Total:

L=Ldata+λphysLphys.\mathcal{L} = \mathcal{L}_{\text{data}} + \lambda_{\text{phys}} \

mathcal{L}_{\text{phys}}.L=Ldata+λphysLphys.

Hyperparameters initial suggestions:

 lambda_phys = 0.1 (start small; tune via validation)

 optimizer: Adam(lr=1e-3) then reduce to 1e-4

 epochs: 100–1000 as needed

 batch size: 64

 weight decay: 1e-5 (optional)

 gradient clipping: 1.0 or lower if instability

4.4 Training loop outline

1. Load training set X, y (y = Teq measured). Fit StandardScaler on X; scale X. Fit scaler on y
(or compute y_mean,y_std) and scale y.

2. Build PINN_Teq(input_dim=D). Move to device cpu or cuda.

3. Optimizer = Adam(model.parameters(), lr).

4. For each epoch:

o iterate mini-batches:

 forward: mu_norm, logvar = model(X_batch) then mu = mu_norm*y_std

+ y_mean, sigma = sqrt(exp(logvar))*y_std

 compute data_nll per sample as above

 compute physics residuals using mu (converted into same units expected

by physics formula) — compute phys_loss = mean(residual^2) on the
batch (Note: compute residuals with tensor operations for autograd; no
detach).

 loss = data_nll.mean() + lambda_phys * phys_loss

 loss.backward(); optional clip; optimizer.step()

o report epoch stats: train data loss, phys loss, validation loss if available.

5. After training, save model.state_dict(), scaler(s), and training history.

Important: phys_loss must be computed on the same samples as data or on a separate
collocation set (recommended to include extra collocation points across P sweep).

5) Inference & plotting

5.1 Inference

 Given a pressure sweep P_sweep (e.g., np.linspace(P_min, P_max, 3000)), build a grid of
input features:

o Use median values for non-pressure features (template medians).

o For the P_bar column, set to the values in the sweep.

 Scale the grid with the saved scaler.

 mu_norm, logvar = model(X_grid_scaled).

 Convert: mu = mu_norm * y_std + y_mean. std = sqrt(exp(logvar))*y_std.

 Produce Pred_Teq = mu, Pred_std = std.

 Save CSV with Pressure, Pred_Teq, Pred_std, Teq_from_curve.

5.2 Plots

PVTSim-style: X-axis Temperature, Y-axis Pressure (many people prefer Teq decreasing
left→right). Steps:

 Plot physical Teq (curve) as plot(Temperature, Pressure) sorted by Temperature (or

Pressure).

 Overlay predicted Teq as plot(Pred_Teq, Pressure) (same sorting), optionally dashed line.

 Add confidence band by fill_betweenx(Pressure, mu - 2*std, mu + 2*std, alpha=0.2).

 Add diagnostics caption: RMSE, MAE, MaxErr, Pearson r.

 Save images hydrate_curve_PVTSim.png, hydrate_curve_overlay.png.

Ensure units are consistent (Temperature °C vs K) before plotting.

6) Diagnostics & tests (what the dev must run & return)

1. Unit checks: print ranges and confirm no columns swapped (Pressure range should be
physically plausible; Temperature range must look like Teq — usually above -40°C and
 less than 300°C). If Temperature max > 200 then it’s probably Kelvin — convert to °C. If
Pressure values are suspicious (e.g., negative where not expected), flag it and attempt
auto-swap if an algorithm detects this (safe corrections only).

2. Tiny forward pass test: run model on 4 sample inputs and confirm shapes: mu.shape ==
(4,1), logvar.shape == (4,1).

3. Physics test: evaluate physics residual on a small P-sample; make sure residual is
numeric and finite. If residual fails (shape mismatch), check broadcasting/reshaping.

4. Full sweep diagnostics: compute RMSE and MAE for the sweep compared to PVTSim Teq
if available. Save metrics.csv.

5. Plot tests: produce both PVTSim-limited view (Temperature X limited to -10..20) and full-
range overlay.

Returnables:

 final_notebook.ipynb (runnable top→bottom)

 pinn_model_teq.pt (or .pth)

 teq_scaler.pkl

 predictions_csv and plotted images

 README.md with usage

7) Uncertainty / Confidence reporting (recommended)

Two approaches:

A — Probabilistic head (preferred)

 Model outputs mu and logvar per sample.

 Interprete std = sqrt(exp(logvar)).

 Use ±2σ band as confidence interval; additionally compute a per-point confidence score
= 1 / (1 + std) or scaled 0..1, high = more confident.

Advantages: single-model, direct NLL loss.

B — Model ensemble
  Train N independent models with different seeds or bootstrap the data.

 Use sample mean as mu, sample std as std.

 Ensemble provides epistemic uncertainty (model uncertainty). Combine with aleatoric

from model (if used).

C — Heuristic risk flag

 If Pred_Teq at pressure lies outside training-range by more than delta (e.g., P outside
training P-range, or Pred_std > threshold), flag as low confidence.

Return in CSV:

7) Uncertainty / Confidence reporting (recommended)

Two approaches:

A — Probabilistic head (preferred)

 Model outputs mu and logvar per sample.

 Interprete std = sqrt(exp(logvar)).

 Use ±2σ band as confidence interval; additionally compute a per-point confidence score
= 1 / (1 + std) or scaled 0..1, high = more confident.

Advantages: single-model, direct NLL loss.

B — Model ensemble

 Train N independent models with different seeds or bootstrap the data.

 Use sample mean as mu, sample std as std.

 Ensemble provides epistemic uncertainty (model uncertainty). Combine with aleatoric

from model (if used).

C — Heuristic risk flag

 If Pred_Teq at pressure lies outside training-range by more than delta (e.g., P outside
training P-range, or Pred_std > threshold), flag as low confidence.

Return in CSV:

Pressure, Pred_Teq, Pred_std, ConfidenceScore (0-1), RiskFlag (0/1)

8) Recommended numeric/engineering defaults (starting point)

 Input features: 8 features (time_s, P_bar, T_K, Sg, phase_frac, pipe_d_m, holdup_liq,
DeltaT)

 Hidden layers: 4 × 128 units

 Activation: ReLU

 Optimizer: Adam lr=1e-3 -> 1e-4

 Batch size: 64

 Epochs: 200–800 (monitor validation)

 λ_phys (physics weight): tune from [1e-3, 1e-2, 1e-1, 1.0]

 Use learning rate scheduler (ReduceLROnPlateau) for training stability

 Save best model by validation RMSE

9) Common pitfalls & how to fix them (from what happened in your notebook)

 Weights-only torch.load issue: When loading PyTorch checkpoints, some saved objects
(pickled objects, older versions) cause UnpicklingError if weights_only semantics
changed. Save a plain state_dict() and load with model.load_state_dict(state_dict,
strict=False). If you must load a full pickled module, use torch.serialization.safe_globals
or re-save from a trusted environment.

 Scaler n_features_in_ mismatch: Make sure the scaler was fit on the same column order
and count used at inference. Save and document X_cols.

 Columns swapped (Pressure vs Temperature): detect by ranges and auto-swap only if

safe. Print a clear warning and save a fixed CSV.

 Units (Kelvin vs °C): detect by max value (>200) and consistently convert.

 Broadcast / shape mismatch when computing residuals: ensure residuals compute with
tensors of compatible shapes; batched operations must match (e.g., use .reshape(-1,1)).

 Extrapolation leads to huge predicted Teq or linear ramp: model may have learned
linear relation or unrealistic extrapolation. Keep training P-range similar to inference P-
range or enforce physics strongly (increase lambda_phys) or augment with synthetic
training points.
  Kernel dying / out-of-memory: reduce batch size, reduce sweep size during diagnostics,
or run on CPU if GPU out-of-memory.

10) Example minimal CSV to ship with handoff (developer test data)

Save this file as example_training_data.csv for dev dry-run (synthetic but sane numbers):

time_s,P_bar,T_K,Sg,phase_frac,pipe_d_m,holdup_liq,DeltaT

0,4.6,293.15,0.198,0.1,0.1,0.05,125.99

6,5.0,293.15,0.2,0.12,0.1,0.06,125.9

12,6.0,293.15,0.197,0.11,0.1,0.05,125.8

...

11) Suggested developer tasks & timeline

1. Day 1: read the paper’s Methodology; implement Teq_from_P via Towler & Mokhatab
OR use PVTSim CSV interpolation. Implement data loader + scaler; validate shapes and
unit conversions.

2. Day 2: implement PINN model and tiny training loop on synthetic dataset; add
probabilistic head; confirm small forward pass and loss decreases on fast toy training.

3. Day 3: add physics residual (implement Bishnoi residual); tune lambda_phys on

validation set; produce predictions vs curve plots and diagnostics (RMSE, MAE).

4. Day 4: add uncertainty reporting (mu,std or ensemble); implement risk flagging; finalize
plotting style (PVTSim).

5. Deliverables: working final_pinn.ipynb, requirements.txt, pinn_model.pt, scaler.pkl,

sample predictions.csv, and README.md.

12) Exact outputs you (the project owner) should expect back

 final_pinn_hydrate.ipynb — runnable top-to-bottom with one cell per step (clean, no

duplicate conflicting cells).

 pinn_model_teq.pt and teq_scaler.pkl (saved artifacts).

 hydrate_curve_prediction_teq.csv containing
Pressure,Teq_curve,Pred_Teq,Pred_std,ConfidenceFlag.
  Plots: hydrate_curve_PVTSim.png, hydrate_curve_overlay.png,
residuals_vs_pressure.png.

 README.md describing how to run (python env, packages), data format, and
parameters.

13) Helpful code snippets (to include in the notebook)

load scaler & model

import pickle, torch

scaler = pickle.load(open('teq_scaler.pkl','rb')) # contains n_features_in_ and mean/std

model = PINN_Teq(input_dim=scaler.n_features_in_)

state = torch.load('pinn_model_teq.pt', map_location='cpu')

model.load_state_dict(state['state_dict'])

model.eval()

build sweep grid & predict

P_sweep = np.linspace(P_min, P_max, 3000)

# build X_grid using medians for non-P features:

templ = {...} # medians for each X_col in same order as scaler

grid = np.array([ row_with_pressure(p, templ) for p in P_sweep ])

X_scaled = scaler.transform(grid)

with torch.no_grad():

mu_norm, logvar = model(torch.tensor(X_scaled, dtype=torch.float32))

mu = mu_norm.numpy().ravel() * y_std + y_mean

std = np.sqrt(np.exp(logvar.numpy().ravel())) * y_std

build sweep grid & predict

P_sweep = np.linspace(P_min, P_max, 3000)

# build X_grid using medians for non-P features:

templ = {...} # medians for each X_col in same order as scaler

grid = np.array([ row_with_pressure(p, templ) for p in P_sweep ])

X_scaled = scaler.transform(grid)

with torch.no_grad():

mu_norm, logvar = model(torch.tensor(X_scaled, dtype=torch.float32))

mu = mu_norm.numpy().ravel() * y_std + y_mean

std = np.sqrt(np.exp(logvar.numpy().ravel())) * y_std

compute residual (example skeleton)

# residual R = f(T_model, P, features) as per Bishnoi-derived expression in paper

# e.g., R = A*(T_model - T_eq_formula(P, features)) -- developer to replace with exact formula

residual = compute_Bishnoi_residual(P_batch, mu_batch, features_batch, constants)

phys_loss = (residual**2).mean()

14) References & where to look in the uploaded PDF

 Methodology / Physical Model: look at Section titled “Methodology” and “Physical

Model for Hydrate Growth Rate” (around pages 5–7). This has the Bishnoi derivation and
references to Towler & Mokhatab and Biberg. The exact algebraic residual expression is
shown there — copy it verbatim into the implementation to compute R.

 Towler & Mokhatab: paper is cited; developer can either copy formula from the paper
(if reproduced) or implement a lookup/interpolation based on PVTSim curve supplied.
15) Final notes / advice for your developer

 Keep everything reproducible: lock exact random seeds, save scalers and model state,
and include requirements.txt (PyTorch version, numpy, pandas, matplotlib, scikit-learn).

 Start small: verify tiny forward passes and physics residual shapes before large sweeps.
Use small sweep sizes for debugging.

 Log extensively: print P/T ranges, n_features_in_, X_cols, and final RMSEs. If column
ordering mismatches, abort and print a red error.

 Unit tests: include a small unit test script that loads sample CSV, runs inference on a 10-
point sweep, and outputs the 10 predictions to a CSV.