Physics Informed Neural Network (Part 2) — Testing the Hypothesis

That’s not a furnace. Image generated by Gemini

Alright. Let’s continue where we left off. Last time we discuss on the theory of Physics Informed Neural Network. We know that it is a method to incorporate physical laws with universal approximation theorem. By putting our first principal equation as loss, we can model certain parts of the equation as trainable parameters. Doing this will help us approximate as close as possible with the actual conditions compare to the best approximation via first principle.

Modify the Loss Function

During my last experiment, we model the loss function based on simple thermodynamic system of a furnace

However, appararently, I did make a mistake. (That 10 years of experience still can’t prevent you from making mistakes apparently). The data that I have collected doesn’t really suitable for the loss function above. The loss function above depends on the temperature of the combustion gas. But typically in a refinery, we don’t have sensor on the combustion gas side. We only have temperature sensor on the inlet and outlet of the fluid heated inside the furnace.

So instead the loss function must be changed. We must focus on the heat absorbed by the process fluid. Following the Laws of Thermodynamics, we assume that the heat produce by the combustion gas will be absorbed by the process fluid leading to the temperature change. But, we must also be realistic. Not all heat released during combustion will be absorbed. Some will be loss. The loss heat must also be taken into account.

This leads to a more accurate differential equation for our system:

Where the terms are now defined based on the available sensor data:

• C(th) is the effective thermal mass of the furnace system (the metal, refractory, etc.), which is a constant we aim to discover.

• Q(released) is the heat generated by fuel gas, modeled as a function of Fuel Gas Flow.

• Q(absorbed) is the heat transferred to the process fluid, which we can calculate directly from our sensors.

• Q(loss) is the heat lost to the surroundings, modeled as being proportional to the outlet temperature.

By replacing the true outlet temperature, T(out) , with the network’s prediction, T(pred), we can formulate a new physics loss. This loss function is a true representation of the system’s thermodynamics as observed by the available sensors, forming a much stronger foundation for our Physics-Informed Neural Network:

Now, we have the loss function that is suitable with the data that we have. For this article, I have already prepared a generic synthetic data that mirror the refinery furnace operation.

You can download the data and the entire notebook from my GitHub repository here:

https://github.com/maercaestro/pinn-experiment

In the dataset we have the

1. Timestep

2. Inlet and Outlet Temperature

3. Inlet Flow

4. Fuel Gas Molecular Weight

5. Excess O2 (excess oxygen)

6. Fuel gas flow

7. Air fuel ratio

To simplify our training we will only be needing the timestep, inlet and outlet temperature, inlet flow (referring to the fluid medium) and fuel gas flow. Excess oxygen and air fuel ratio will contribute to the furnace efficiency, something that we may not need now but crucial actually to build a more accurate model.

Let’s load the data, and prepare it for our training. We will load it via pandas, normalize it using MinMaxScaler and load it as tensor.

#first, let's import the libraries and packages
from sklearn.preprocessing import MinMaxScaler
import torch
import pandas as pd

#load the dataset
combined_df = pd.read_csv('curated_training_data.csv')
combined_df.info()

If the data is loaded correctly, you will see results as below

The data has been specially curated out of combine steady-state and transient condition of a furnace. Therefore, all the major preprocessing has been done. Aas you can see above we already have the timestep and all the required features. But, we did have extra date columns, let’s remove that.

features_to_scale = combined_df.drop(columns=['Date'])

#now it's time to scaler the data
scaler = MinMaxScaler()

# Fit the scaler on the data and transform it
data_scaled = scaler.fit_transform(features_to_scale)

print("Data has been successfully scaled.")
print("Shape of scaled data:", data_scaled.shape)

The data has been scaled. Now, let’s convert our data to tensor and get it ready for training.

# we need to convert our data to tensor
data_tensor = torch.tensor(data_scaled, dtype=torch.float32)

# Get the column names from the scaled dataframe for easy indexing
scaled_columns = features_to_scale.columns

# The input to the NN is time 't'
t = data_tensor[:, scaled_columns.get_loc('t')].view(-1, 1)

# The target for the data loss is 'OutletTemp'
T_outlet = data_tensor[:, scaled_columns.get_loc('OutletTemp')].view(-1, 1)

# A dictionary to hold all the features needed for the physics loss calculation
physics_features = {
    'InletTemp': data_tensor[:, scaled_columns.get_loc('InletTemp')].view(-1, 1),
    'InletFlow': data_tensor[:, scaled_columns.get_loc('InletFlow')].view(-1, 1),
    'FuelGasFlow': data_tensor[:, scaled_columns.get_loc('FuelGasFlow')].view(-1, 1),
    'AirFuelRatio': data_tensor[:, scaled_columns.get_loc('AirFuelRatio')].view(-1, 1),
    'ExcessO2': data_tensor[:, scaled_columns.get_loc('ExcessO2')].view(-1, 1),
    'FuelGasMW': data_tensor[:, scaled_columns.get_loc('FuelGasMW')].view(-1, 1),
}

print("\nData successfully converted to PyTorch tensors.")
print("Shape of time tensor 't':", t.shape)
print("Shape of target tensor 'T_outlet':", T_outlet.shape)
print("Feature 'InletTemp' tensor shape:", physics_features['InletTemp'].shape)

If it is succesful, you will get an output like below:

Now, our data is ready, let’s build our model and our loss functions. We will build a simple network with 3 linear layers, followed by Tanh activation functions. Inside the network, we will add our learnable parameters as log of the parameters to ensure the model can effectively learn the reduce the gap between thermodynamic properties and actual condition. Originally, I try training this without the log of the parameters. However, it proves troublesome as it introduce negative values to the parameters leading to unstable network.

So we change our strategy by defining our parameters as log of the physical constants. Then we calculate the coefficients by taking the explonents. This gives mathematical guarantee that our network will be stable during training.

import torch
import torch.nn as nn
import matplotlib.pyplot as plt
import pandas as pd

# --- 1. Improved Model Definition (with Log-Parameters) ---
class FurnacePINN(nn.Module):
    def __init__(self):
        super(FurnacePINN, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(1, 64), nn.Tanh(),
            nn.Linear(64, 64), nn.Tanh(),
            nn.Linear(64, 1)
        )
        # We learn the log of the parameters, initializing at 0.0 (so the param starts at exp(0)=1.0)
        self.log_C_th = nn.Parameter(torch.tensor([0.0]))
        self.log_Cp_naphtha = nn.Parameter(torch.tensor([0.0]))
        self.log_k_LHV = nn.Parameter(torch.tensor([0.0]))
        self.log_k_loss = nn.Parameter(torch.tensor([0.0]))
        
        min_temp = scaler.data_min_[scaled_columns.get_loc('OutletTemp')]
        max_temp = scaler.data_max_[scaled_columns.get_loc('OutletTemp')]
        scaled_T_ambient = (298.15 - min_temp) / (max_temp - min_temp)
        self.T_ambient = torch.tensor(scaled_T_ambient, dtype=torch.float32)

    def forward(self, t):
        return self.net(t)

Now let’s build the loss function

#let's build the loss function
def compute_loss(model, t, T_outlet, physics_features):
    T_pred = model(t)
    loss_data = torch.mean((T_pred - T_outlet)**2)
    
    t.requires_grad = True
    T_pred_phys = model(t)
    dT_dt = torch.autograd.grad(T_pred_phys, t, grad_outputs=torch.ones_like(T_pred_phys), create_graph=True)[0]
    
    # Exponentiate the logs to get guaranteed-positive physical parameters
    C_th = torch.exp(model.log_C_th)
    Cp_naphtha = torch.exp(model.log_Cp_naphtha)
    k_LHV = torch.exp(model.log_k_LHV)
    k_loss = torch.exp(model.log_k_loss)
    
    T_in, m_in, F_in = physics_features['InletTemp'], physics_features['InletFlow'], physics_features['FuelGasFlow']
    
    Q_released = k_LHV * F_in
    Q_absorbed = m_in * Cp_naphtha * (T_pred_phys - T_in)
    Q_loss = k_loss * (T_pred_phys - model.T_ambient)
    
    residual = C_th * dT_dt - (Q_released - Q_absorbed - Q_loss)
    loss_phys = torch.mean(residual**2)
    
    total_loss = loss_data + 1.0 * loss_phys
    return total_loss, loss_data, loss_phys

Once it’s done, let’s define our training loop. As usual, we will use Adam as optimizer, with learning rate of 1e-3. We will train this at 2000 epochs. Don’t forget to log the training loss so we can use it to plot our data.

model = FurnacePINN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
epochs = 20000

# Lists to store log data
log_epochs = []
log_data_loss = []
log_phys_loss = []
log_params = {'C_th': [], 'Cp_n': [], 'k_LHV': [], 'k_loss': []}

print("\n--- Starting Model Training ---")
for epoch in range(epochs + 1):
    optimizer.zero_grad()
    loss, loss_d, loss_p = compute_loss(model, t, T_outlet, physics_features)
    loss.backward()
    optimizer.step()
    
    # Logging data at specified intervals
    if epoch % 100 == 0:
        log_epochs.append(epoch)
        log_data_loss.append(loss_d.item())
        log_phys_loss.append(loss_p.item())
        
        # Log the actual (exponentiated) parameter values
        log_params['C_th'].append(torch.exp(model.log_C_th).item())
        log_params['Cp_n'].append(torch.exp(model.log_Cp_naphtha).item())
        log_params['k_LHV'].append(torch.exp(model.log_k_LHV).item())
        log_params['k_loss'].append(torch.exp(model.log_k_loss).item())
        
        if epoch % 1000 == 0:
            print(f"Epoch {epoch}: Loss={loss.item():.6f}")

print("\n--- Training Complete ---")

Once the training completed, we can use the log data to plot our training loops.

# --- 4. Immediate Visualization ---
print("--- Generating Plots ---")
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 10), sharex=True)
fig.suptitle('Training Progress Analysis (with Positive Constraints)', fontsize=16)

# Plot 1: Loss Curves
ax1.plot(log_epochs, log_data_loss, label='Data Loss', color='blue')
ax1.plot(log_epochs, log_phys_loss, label='Physics Loss', color='red')
ax1.set_yscale('log')
ax1.set_ylabel('Loss (Log Scale)')
ax1.set_title('Loss Curves')
ax1.legend()
ax1.grid(True, which="both", ls="--")

# Plot 2: Parameter Trajectories
ax2.plot(log_epochs, log_params['C_th'], label='C_th (Thermal Mass)')
ax2.plot(log_epochs, log_params['Cp_n'], label='Cp_n (Specific Heat)')
ax2.plot(log_epochs, log_params['k_LHV'], label='k_LHV (Heating Value)')
ax2.plot(log_epochs, log_params['k_loss'], label='k_loss (Heat Loss)')
ax2.set_ylabel('Parameter Value')
ax2.set_xlabel('Epoch')
ax2.set_title('Parameter Trajectories')
ax2.legend()
ax2.grid(True, ls="--")

fig.tight_layout(rect=[0, 0, 1, 0.96])
plt.savefig('final_training_analysis.png')
print("Generated final plot: 'final_training_analysis.png'")

This will output our graph as below:

We see that both the data loss and physics loss reduce over time. And that is a good sign. And we also achieve a good convergance of parameter value when our data loss and physics loss reach low level.

Now, it’s time for the true test. How well our model behave on their prediction. Does it really learn? So let’s put our predicted value on top of the actual curve

import numpy as np
import matplotlib.pyplot as plt
import torch

# --- 1. Get Model Predictions ---
# Set the model to evaluation mode (disables things like dropout, not used here but good practice)
model.eval() 

# Disable gradient calculations for faster inference
with torch.no_grad(): 
    T_pred_scaled = model(t)

print("--- Step 1: Model predictions generated. ---")

# --- 2. Inverse Transform Data to Real Units ---
# Convert the PyTorch tensors to NumPy arrays for use with scikit-learn and Matplotlib
T_actual_scaled_np = T_outlet.numpy()
T_pred_scaled_np = T_pred_scaled.numpy()
# Use the unscaled time from the original combined_df for a more readable x-axis
time_real_units = combined_df['t'].to_numpy() 

# The scaler was fitted on a DataFrame with multiple columns. To inverse transform,
# we need to create a dummy array of the same shape and put our temperature data
# in the correct column.
num_features = len(scaler.data_min_)
outlet_temp_col_index = list(scaled_columns).index('OutletTemp')

# Create a dummy array for the actual temperature values
actual_full_shape = np.zeros((len(T_actual_scaled_np), num_features))
actual_full_shape[:, outlet_temp_col_index] = T_actual_scaled_np.flatten()

# Create a dummy array for the predicted temperature values
pred_full_shape = np.zeros((len(T_pred_scaled_np), num_features))
pred_full_shape[:, outlet_temp_col_index] = T_pred_scaled_np.flatten()

# Use the scaler to inverse transform the data back to its original scale
T_actual_real = scaler.inverse_transform(actual_full_shape)[:, outlet_temp_col_index]
T_pred_real = scaler.inverse_transform(pred_full_shape)[:, outlet_temp_col_index]

print("--- Step 2: Predictions converted back to real units. ---")

# --- 3. Plot the Final Results ---
plt.figure(figsize=(15, 7))
plt.plot(time_real_units, T_actual_real, label='Actual Temperature', color='black', linewidth=2.5)
plt.plot(time_real_units, T_pred_real, label='PINN Predicted Temperature', color='red', linestyle='--', linewidth=2)

plt.xlabel('Time (s)')
plt.ylabel('Outlet Temperature (degC)')
plt.title('Final PINN Performance Evaluation', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, which="both", ls="--", alpha=0.6)
plt.show()

print("\n--- Step 3: Visualization complete. ---")

# --- 4. Display Final Learned Parameters ---
print("\nFinal Learned Physical Parameter Coefficients:")
print(f"  - C_th (Thermal Mass): {torch.exp(model.log_C_th).item():.4f}")
print(f"  - Cp_n (Specific Heat): {torch.exp(model.log_Cp_naphtha).item():.4f}")
print(f"  - k_LHV (Heating Value): {torch.exp(model.log_k_LHV).item():.4f}")
print(f"  - k_loss (Heat Loss): {torch.exp(model.log_k_loss).item():.4f}")

If everything works correctly, we will get a good graph like this

Which is actually a good results, for a first try. Bear in mind, that we haven’t fully split our training data. So an overfitting case might happen here. But I welcome you try that part yourself, and maybe tell me what happen when we split between training and testing.

For now, this is PINN, a neural network that understand physics. At the end of the code, it will output the learnable parameters. That hasn’t been scaled actually. So, you have to scale that back.

That’s all for now. Thanks for reading, and see you next time.