Developer Handoff — PINN for Hydrate Equilibrium (Teq) prediction
Source: OTC/SPE paper supplied ([Link]). 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)
[Link] 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([Link]):
def __init__(self, input_dim):
super().__init__()
[Link] = [Link](
[Link](input_dim, 128), [Link](),
[Link](128, 128), [Link](),
[Link](128, 128), [Link](),
[Link](128, 128), [Link](),
self.mu_head = [Link](128, 1)
self.logvar_head = [Link](128, 1)
def forward(self, x):
h = [Link](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([Link](), 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
[Link](); optional clip; [Link]()
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., [Link](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: [Link] ==
(4,1), [Link] == (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 [Link].
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
[Link] 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 [Link] 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 [Link].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, [Link], pinn_model.pt, [Link],
sample [Link], and [Link].
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.
[Link] 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 = [Link](open('teq_scaler.pkl','rb')) # contains n_features_in_ and mean/std
model = PINN_Teq(input_dim=scaler.n_features_in_)
state = [Link]('pinn_model_teq.pt', map_location='cpu')
model.load_state_dict(state['state_dict'])
[Link]()
build sweep grid & predict
P_sweep = [Link](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 = [Link]([ row_with_pressure(p, templ) for p in P_sweep ])
X_scaled = [Link](grid)
with torch.no_grad():
mu_norm, logvar = model([Link](X_scaled, dtype=torch.float32))
mu = mu_norm.numpy().ravel() * y_std + y_mean
std = [Link]([Link]([Link]().ravel())) * y_std
build sweep grid & predict
P_sweep = [Link](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 = [Link]([ row_with_pressure(p, templ) for p in P_sweep ])
X_scaled = [Link](grid)
with torch.no_grad():
mu_norm, logvar = model([Link](X_scaled, dtype=torch.float32))
mu = mu_norm.numpy().ravel() * y_std + y_mean
std = [Link]([Link]([Link]().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 [Link] (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.