Heat Transfer#
Heat Transfer Example for JAXSR.
Demonstrates discovering heat transfer correlations from data, including:
Nusselt number correlations
Natural convection
Forced convection
import jax.numpy as jnp
import numpy as np
from jaxsr import BasisLibrary, Constraints, SymbolicRegressor
Discover Dittus-Boelter correlation for turbulent forced convection.#
True model (Dittus–Boelter):
\[\mathrm{Nu} = 0.023\,\mathrm{Re}^{0.8}\,\mathrm{Pr}^{0.4}\]
np.random.seed(42)
n_samples = 100
# Reynolds and Prandtl number ranges for turbulent flow
Re = np.random.uniform(10000, 100000, n_samples)
Pr = np.random.uniform(0.7, 100, n_samples)
# Dittus-Boelter correlation
Nu_true = 0.023 * Re**0.8 * Pr**0.4
Nu = Nu_true * (1 + np.random.randn(n_samples) * 0.05)
# Work in log space for power law discovery
log_Re = np.log(Re)
log_Pr = np.log(Pr)
log_Nu = np.log(Nu)
X = jnp.column_stack([log_Re, log_Pr])
y = jnp.array(log_Nu)
print("\nTrue model: Nu = 0.023 * Re^0.8 * Pr^0.4")
print(f"Log form: ln(Nu) = {np.log(0.023):.3f} + 0.8*ln(Re) + 0.4*ln(Pr)")
library = (
BasisLibrary(n_features=2, feature_names=["ln_Re", "ln_Pr"])
.add_constant()
.add_linear()
.add_interactions(max_order=2)
)
model = SymbolicRegressor(
basis_library=library,
max_terms=4,
strategy="greedy_forward",
)
model.fit(X, y)
print("\nDiscovered expression (log space):")
print(f" {model.expression_}")
print(f" R² = {model.metrics_['r2']:.4f}")
# Extract exponents
if "ln_Re" in model.selected_features_:
idx = model.selected_features_.index("ln_Re")
re_exp = float(model.coefficients_[idx])
print(f"\nRe exponent: {re_exp:.2f} (true: 0.80)")
if "ln_Pr" in model.selected_features_:
idx = model.selected_features_.index("ln_Pr")
pr_exp = float(model.coefficients_[idx])
print(f"Pr exponent: {pr_exp:.2f} (true: 0.40)")
print("\n--- Why the exponents differ from true values ---")
print("The spurious ln_Re*ln_Pr interaction term absorbs variance that should")
print("belong to individual terms, biasing the recovered exponents.")
if "ln_Re" in model.selected_features_:
idx = model.selected_features_.index("ln_Re")
re_exp = float(model.coefficients_[idx])
print(f" - Re: {re_exp:.2f} vs true 0.80 (error: {abs(re_exp - 0.8)/0.8*100:.1f}%)")
if "ln_Pr" in model.selected_features_:
idx = model.selected_features_.index("ln_Pr")
pr_exp = float(model.coefficients_[idx])
print(f" - Pr: {pr_exp:.2f} vs true 0.40 (error: {abs(pr_exp - 0.4)/0.4*100:.1f}%)")
print()
print("Strategies for exact recovery:")
print(" 1. Exclude interaction terms from basis library (if power law is known)")
print(" 2. Use exhaustive search instead of greedy forward selection")
print(" 3. Use parametric basis functions: add_parametric() with profile likelihood")
print(" (see langmuir_doe_active_learning.ipynb for example)")
# Parameter significance, diagnostics, and ANOVA
import matplotlib.pyplot as plt
from scipy import stats as sp_stats
from jaxsr import anova
from jaxsr.plotting import plot_parity
intervals = model.coefficient_intervals(alpha=0.05)
n_obs, k_terms = len(np.asarray(y)), len(model.selected_features_)
df_resid = n_obs - k_terms
print("Parameter Significance (95% CI):")
print(f" {'Term':>15s} {'Estimate':>10s} {'Std Err':>9s} {'t':>8s} {'p-value':>10s} 95% CI")
print(" " + "-" * 75)
for name, (est, lo, hi, se) in intervals.items():
t_val = est / se if abs(se) > 1e-15 else float("inf")
p_val = float(2 * (1 - sp_stats.t.cdf(abs(t_val), df_resid))) if df_resid > 0 else 0.0
sig = "***" if p_val < 0.001 else ("**" if p_val < 0.01 else ("*" if p_val < 0.05 else ""))
print(f" {name:>15s} {est:10.4f} {se:9.4f} {t_val:8.2f} {p_val:10.2e} [{lo:.4f}, {hi:.4f}] {sig}")
print(" --- *** p<0.001, ** p<0.01, * p<0.05")
# Parity and residual plots
y_pred = model.predict(X)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_parity(y, y_pred, ax=axes[0], title="Forced Convection: Parity")
residuals = np.array(y - y_pred)
axes[1].scatter(np.array(y_pred), residuals, alpha=0.6, c="steelblue", edgecolors="white", linewidth=0.5)
axes[1].axhline(y=0, color="r", linestyle="--")
axes[1].set_xlabel("Predicted")
axes[1].set_ylabel("Residuals")
axes[1].set_title("Forced Convection: Residuals")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# ANOVA
anova_result = anova(model)
summary_sources = {"Model", "Residual", "Total"}
print("\nANOVA Table (Forced Convection)")
print("=" * 80)
print(f" {'Source':25s} {'DF':>4} {'Sum Sq':>12} {'Mean Sq':>12} {'F':>10} {'p-value':>10}")
print("-" * 80)
for row in anova_result.rows:
f_str = f"{row.f_value:10.2f}" if row.f_value is not None else " "
p_str = f"{row.p_value:10.4f}" if row.p_value is not None else " "
print(f" {row.source:25s} {row.df:4d} {row.sum_sq:12.4f} {row.mean_sq:12.4f} {f_str} {p_str}")
print("-" * 80)
term_rows = [r for r in anova_result.rows if r.source not in summary_sources]
if term_rows:
model_ss = sum(r.sum_sq for r in term_rows)
print("\nVariance Contributions:")
print("(Percentages relative to Model SS, not Total SS — shows relative importance within fitted model)")
for row in term_rows:
pct = 100 * row.sum_sq / model_ss if model_ss > 0 else 0
sig = (
"***" if row.p_value is not None and row.p_value < 0.001 else (
"**" if row.p_value is not None and row.p_value < 0.01 else (
"*" if row.p_value is not None and row.p_value < 0.05 else ""))
)
print(f" {row.source:25s} {pct:6.1f}% {sig}")
Interpretation: Dittus-Boelter exponent recovery#
The discovered expression recovers the power-law structure of the Dittus-Boelter correlation but with noticeable deviations in the exponents:
Re exponent: 0.77 vs. true 0.80 (approximately 4% low)
Pr exponent: 0.34 vs. true 0.40 (approximately 15% low)
The model also includes a spurious ln_Re*ln_Pr interaction term. In the original
(non-log) space, this interaction corresponds to a Re^a * Pr^b coupling that is
not present in the true correlation. This term absorbs some of the variance that
should be attributed to the individual ln_Re and ln_Pr terms, pulling their
coefficients (i.e., the recovered exponents) away from the true values.
This is a known behavior of greedy forward selection: the interaction term may be selected early because it greedily reduces the residual sum of squares before the correct individual terms are fully established at their true coefficients. Once the interaction term is in the model, the remaining terms adjust their coefficients to compensate, leading to biased exponent estimates.
Strategies for exact recovery:
Use a more constrained basis library containing only
ln_Re,ln_Pr, and a constant (excluding interaction terms), since the true model is known to be a simple power law.Use exhaustive search instead of greedy forward selection to evaluate all candidate subsets.
Use parametric (nonlinear) basis functions to fit
Re^aandPr^bdirectly, rather than relying on the log-space linearization.
Discover natural convection correlation.#
True model (vertical plate, laminar):
\[\mathrm{Nu} = 0.59\,\mathrm{Ra}^{0.25}\]
np.random.seed(42)
n_samples = 80
# Rayleigh number range for laminar flow
Ra = np.random.uniform(1e4, 1e9, n_samples)
# Churchill-Chu correlation (simplified for laminar)
C = 0.59
n = 0.25
Nu_true = C * Ra**n
Nu = Nu_true * (1 + np.random.randn(n_samples) * 0.03)
# Log transformation
log_Ra = np.log10(Ra)
log_Nu = np.log10(Nu)
X = jnp.array(log_Ra).reshape(-1, 1)
y = jnp.array(log_Nu)
print("\nTrue model: Nu = 0.59 * Ra^0.25")
print(f"Log form: log10(Nu) = {np.log10(C):.3f} + 0.25*log10(Ra)")
library = BasisLibrary(n_features=1, feature_names=["log_Ra"]).add_constant().add_linear()
model = SymbolicRegressor(
basis_library=library,
max_terms=2,
strategy="exhaustive",
)
model.fit(X, y)
print("\nDiscovered expression:")
print(f" {model.expression_}")
print(f" R² = {model.metrics_['r2']:.4f}")
True model: Nu = 0.59 * Ra^0.25
Log form: log10(Nu) = -0.229 + 0.25*log10(Ra)
Discovered expression:
y = - 0.2709 + 0.2549*log_Ra
R² = 0.9884
# Diagnostics and ANOVA for natural convection
intervals = model.coefficient_intervals(alpha=0.05)
n_obs, k_terms = len(np.asarray(y)), len(model.selected_features_)
df_resid = n_obs - k_terms
print("Parameter Significance (95% CI):")
print(f" {'Term':>15s} {'Estimate':>10s} {'Std Err':>9s} {'t':>8s} {'p-value':>10s} 95% CI")
print(" " + "-" * 75)
for name, (est, lo, hi, se) in intervals.items():
t_val = est / se if abs(se) > 1e-15 else float("inf")
p_val = float(2 * (1 - sp_stats.t.cdf(abs(t_val), df_resid))) if df_resid > 0 else 0.0
sig = "***" if p_val < 0.001 else ("**" if p_val < 0.01 else ("*" if p_val < 0.05 else ""))
print(f" {name:>15s} {est:10.4f} {se:9.4f} {t_val:8.2f} {p_val:10.2e} [{lo:.4f}, {hi:.4f}] {sig}")
print(" --- *** p<0.001, ** p<0.01, * p<0.05")
y_pred = model.predict(X)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_parity(y, y_pred, ax=axes[0], title="Natural Convection: Parity")
residuals = np.array(y - y_pred)
axes[1].scatter(np.array(y_pred), residuals, alpha=0.6, c="steelblue", edgecolors="white", linewidth=0.5)
axes[1].axhline(y=0, color="r", linestyle="--")
axes[1].set_xlabel("Predicted")
axes[1].set_ylabel("Residuals")
axes[1].set_title("Natural Convection: Residuals")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
anova_result = anova(model)
summary_sources = {"Model", "Residual", "Total"}
print("\nANOVA Table (Natural Convection)")
print("=" * 80)
print(f" {'Source':25s} {'DF':>4} {'Sum Sq':>12} {'Mean Sq':>12} {'F':>10} {'p-value':>10}")
print("-" * 80)
for row in anova_result.rows:
f_str = f"{row.f_value:10.2f}" if row.f_value is not None else " "
p_str = f"{row.p_value:10.4f}" if row.p_value is not None else " "
print(f" {row.source:25s} {row.df:4d} {row.sum_sq:12.4f} {row.mean_sq:12.4f} {f_str} {p_str}")
print("-" * 80)
term_rows = [r for r in anova_result.rows if r.source not in summary_sources]
if term_rows:
model_ss = sum(r.sum_sq for r in term_rows)
print("\nVariance Contributions:")
print("(Percentages relative to Model SS, not Total SS — shows relative importance within fitted model)")
for row in term_rows:
pct = 100 * row.sum_sq / model_ss if model_ss > 0 else 0
sig = (
"***" if row.p_value is not None and row.p_value < 0.001 else (
"**" if row.p_value is not None and row.p_value < 0.01 else (
"*" if row.p_value is not None and row.p_value < 0.05 else ""))
)
print(f" {row.source:25s} {pct:6.1f}% {sig}")
Discover fin efficiency correlation.#
True model:
\[\eta = \frac{\tanh(mL)}{mL}, \quad mL = \sqrt{\frac{hP}{kA}}\,L\]
np.random.seed(42)
n_samples = 100
# mL parameter (dimensionless fin parameter)
mL = np.random.uniform(0.1, 3.0, n_samples)
# True fin efficiency
eta_true = np.tanh(mL) / mL
eta = eta_true + np.random.randn(n_samples) * 0.01
X = jnp.array(mL).reshape(-1, 1)
y = jnp.array(eta)
print("\nTrue model: eta = tanh(mL) / mL")
# Build library with hyperbolic functions
library = (
BasisLibrary(n_features=1, feature_names=["mL"])
.add_constant()
.add_linear()
.add_polynomials(max_degree=4)
.add_transcendental(["tanh", "inv"])
)
# Add the exact form as a custom function
library.add_custom(
name="tanh(mL)/mL",
func=lambda X: jnp.tanh(X[:, 0]) / X[:, 0],
complexity=3,
)
# Constraint: efficiency must be between 0 and 1
constraints = Constraints().add_bounds("y", lower=0.0, upper=1.0)
model = SymbolicRegressor(
basis_library=library,
max_terms=3,
strategy="greedy_forward",
constraints=constraints,
)
model.fit(X, y)
print("\nDiscovered expression:")
print(f" {model.expression_}")
print(f" R² = {model.metrics_['r2']:.4f}")
True model: eta = tanh(mL) / mL
Discovered expression:
y = 1.001*tanh(mL)/mL
R² = 0.9984
# Diagnostics and ANOVA for fin efficiency
intervals = model.coefficient_intervals(alpha=0.05)
n_obs, k_terms = len(np.asarray(y)), len(model.selected_features_)
df_resid = n_obs - k_terms
print("Parameter Significance (95% CI):")
print(f" {'Term':>15s} {'Estimate':>10s} {'Std Err':>9s} {'t':>8s} {'p-value':>10s} 95% CI")
print(" " + "-" * 75)
for name, (est, lo, hi, se) in intervals.items():
t_val = est / se if abs(se) > 1e-15 else float("inf")
p_val = float(2 * (1 - sp_stats.t.cdf(abs(t_val), df_resid))) if df_resid > 0 else 0.0
sig = "***" if p_val < 0.001 else ("**" if p_val < 0.01 else ("*" if p_val < 0.05 else ""))
print(f" {name:>15s} {est:10.4f} {se:9.4f} {t_val:8.2f} {p_val:10.2e} [{lo:.4f}, {hi:.4f}] {sig}")
print(" --- *** p<0.001, ** p<0.01, * p<0.05")
y_pred = model.predict(X)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_parity(y, y_pred, ax=axes[0], title="Fin Efficiency: Parity")
residuals = np.array(y - y_pred)
axes[1].scatter(np.array(y_pred), residuals, alpha=0.6, c="steelblue", edgecolors="white", linewidth=0.5)
axes[1].axhline(y=0, color="r", linestyle="--")
axes[1].set_xlabel("Predicted")
axes[1].set_ylabel("Residuals")
axes[1].set_title("Fin Efficiency: Residuals")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
anova_result = anova(model)
summary_sources = {"Model", "Residual", "Total"}
print("\nANOVA Table (Fin Efficiency)")
print("=" * 80)
print(f" {'Source':25s} {'DF':>4} {'Sum Sq':>12} {'Mean Sq':>12} {'F':>10} {'p-value':>10}")
print("-" * 80)
for row in anova_result.rows:
f_str = f"{row.f_value:10.2f}" if row.f_value is not None else " "
p_str = f"{row.p_value:10.4f}" if row.p_value is not None else " "
print(f" {row.source:25s} {row.df:4d} {row.sum_sq:12.4f} {row.mean_sq:12.4f} {f_str} {p_str}")
print("-" * 80)
term_rows = [r for r in anova_result.rows if r.source not in summary_sources]
if term_rows:
model_ss = sum(r.sum_sq for r in term_rows)
print("\nVariance Contributions:")
print("(Percentages relative to Model SS, not Total SS — shows relative importance within fitted model)")
for row in term_rows:
pct = 100 * row.sum_sq / model_ss if model_ss > 0 else 0
sig = (
"***" if row.p_value is not None and row.p_value < 0.001 else (
"**" if row.p_value is not None and row.p_value < 0.01 else (
"*" if row.p_value is not None and row.p_value < 0.05 else ""))
)
print(f" {row.source:25s} {pct:6.1f}% {sig}")
Discover heat exchanger effectiveness-NTU relationship.#
True model (parallel flow):
\[\varepsilon = \frac{1 - e^{-\mathrm{NTU}(1+C)}}{1+C}\]
np.random.seed(42)
n_samples = 120
# NTU and capacity ratio
NTU = np.random.uniform(0.1, 5.0, n_samples)
C = np.random.uniform(0.0, 1.0, n_samples)
# Parallel flow effectiveness
eps_true = (1 - np.exp(-NTU * (1 + C))) / (1 + C)
eps = eps_true + np.random.randn(n_samples) * 0.01
X = jnp.column_stack([NTU, C])
y = jnp.array(eps)
print("\nTrue model: eps = (1 - exp(-NTU*(1+C))) / (1+C)")
library = (
BasisLibrary(n_features=2, feature_names=["NTU", "C"])
.add_constant()
.add_linear()
.add_polynomials(max_degree=2)
.add_interactions(max_order=2)
.add_transcendental(["exp"])
)
# Add specific forms
library.add_custom(
name="exp(-NTU)",
func=lambda X: jnp.exp(-X[:, 0]),
complexity=2,
)
library.add_custom(
name="exp(-NTU*(1+C))",
func=lambda X: jnp.exp(-X[:, 0] * (1 + X[:, 1])),
complexity=3,
)
library.add_custom(
name="1/(1+C)",
func=lambda X: 1 / (1 + X[:, 1]),
complexity=2,
)
# Effectiveness between 0 and 1
constraints = Constraints().add_bounds("y", lower=0.0, upper=1.0)
model = SymbolicRegressor(
basis_library=library,
max_terms=5,
strategy="greedy_forward",
constraints=constraints,
)
model.fit(X, y)
print("\nDiscovered expression:")
print(f" {model.expression_}")
print(f" R² = {model.metrics_['r2']:.4f}")
print("\n--- Structural Limitation of Linear Basis Approach ---")
print("True form: ε = [1 - exp(-NTU*(1+C))] / (1+C)")
print()
print("This is a PRODUCT of two nonlinear terms, not a linear combination.")
print("A linear SR model can only assign CONSTANT coefficients to each basis function.")
print("True model requires coefficient on exp() term to vary with C as -1/(1+C),")
print("which is not representable in linear-in-parameters form.")
print()
print("Despite this limitation, the polynomial approximation achieves excellent fit")
print("(R² > 0.99) and is practical for interpolation within the sampled range.")
print("Caution: Extrapolation beyond training range may be inaccurate.")
# Diagnostics and ANOVA for heat exchanger
intervals = model.coefficient_intervals(alpha=0.05)
n_obs, k_terms = len(np.asarray(y)), len(model.selected_features_)
df_resid = n_obs - k_terms
print("Parameter Significance (95% CI):")
print(f" {'Term':>15s} {'Estimate':>10s} {'Std Err':>9s} {'t':>8s} {'p-value':>10s} 95% CI")
print(" " + "-" * 75)
for name, (est, lo, hi, se) in intervals.items():
t_val = est / se if abs(se) > 1e-15 else float("inf")
p_val = float(2 * (1 - sp_stats.t.cdf(abs(t_val), df_resid))) if df_resid > 0 else 0.0
sig = "***" if p_val < 0.001 else ("**" if p_val < 0.01 else ("*" if p_val < 0.05 else ""))
print(f" {name:>15s} {est:10.4f} {se:9.4f} {t_val:8.2f} {p_val:10.2e} [{lo:.4f}, {hi:.4f}] {sig}")
print(" --- *** p<0.001, ** p<0.01, * p<0.05")
y_pred = model.predict(X)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
plot_parity(y, y_pred, ax=axes[0], title="Heat Exchanger: Parity")
residuals = np.array(y - y_pred)
axes[1].scatter(np.array(y_pred), residuals, alpha=0.6, c="steelblue", edgecolors="white", linewidth=0.5)
axes[1].axhline(y=0, color="r", linestyle="--")
axes[1].set_xlabel("Predicted")
axes[1].set_ylabel("Residuals")
axes[1].set_title("Heat Exchanger: Residuals")
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
anova_result = anova(model)
summary_sources = {"Model", "Residual", "Total"}
print("\nANOVA Table (Heat Exchanger)")
print("=" * 80)
print(f" {'Source':25s} {'DF':>4} {'Sum Sq':>12} {'Mean Sq':>12} {'F':>10} {'p-value':>10}")
print("-" * 80)
for row in anova_result.rows:
f_str = f"{row.f_value:10.2f}" if row.f_value is not None else " "
p_str = f"{row.p_value:10.4f}" if row.p_value is not None else " "
print(f" {row.source:25s} {row.df:4d} {row.sum_sq:12.4f} {row.mean_sq:12.4f} {f_str} {p_str}")
print("-" * 80)
term_rows = [r for r in anova_result.rows if r.source not in summary_sources]
if term_rows:
model_ss = sum(r.sum_sq for r in term_rows)
print("\nVariance Contributions:")
print("(Percentages relative to Model SS, not Total SS — shows relative importance within fitted model)")
for row in term_rows:
pct = 100 * row.sum_sq / model_ss if model_ss > 0 else 0
sig = (
"***" if row.p_value is not None and row.p_value < 0.001 else (
"**" if row.p_value is not None and row.p_value < 0.01 else (
"*" if row.p_value is not None and row.p_value < 0.05 else ""))
)
print(f" {row.source:25s} {pct:6.1f}% {sig}")
Interpretation: effectiveness-NTU structural limitations#
The discovered model achieves an excellent R-squared of 0.991 but does not recover the exact analytical form. The key discrepancies are:
The coefficient on
1/(1+C)is 0.878 (true: 1.0).The coefficient on
exp(-NTU*(1+C))is a constant -0.627, whereas in the true expression(1 - exp(-NTU*(1+C))) / (1+C), the effective coefficient on the exponential term is-1/(1+C), which varies with C.The model compensates with three additional correction terms (
C,NTU*C,NTU) that approximate the residual error from the structural mismatch.
This is a structural limitation of the linear-combination basis approach. The true form is a product of two nonlinear terms:
\[\varepsilon = \frac{1}{1+C} \cdot \bigl(1 - e^{-\mathrm{NTU}(1+C)}\bigr)\]
This product cannot be exactly represented as a linear combination of the
library’s basis functions, because the coefficient on exp(-NTU*(1+C)) in the
true model is itself a function of C (specifically, -1/(1+C)), not a constant.
A linear symbolic regression model can only assign constant coefficients to each
basis function, so it cannot capture this coupling exactly.
Despite this structural limitation, the R-squared of 0.991 indicates that the
five-term approximation is practically excellent across the sampled NTU and C
ranges. For exact recovery, one would need to include the composite basis
function exp(-NTU*(1+C))/(1+C) directly in the library, or use a nonlinear
regression approach that allows variable-dependent coefficients.