Error on paramater estimation using MCMC

[PapaKwansa]

In this simulation, I am trying to match the simulated dataset, which was generated in COMSOL, with a synthetic dataset in the codes below by using MCMC methods for the posterior parameter estimation. The MCMC algorithm has not been able to successfully access the datasets from my COMSOL model. When I run model.exports() it returns null. This is my first time doing this, and I would appreciate your help and suggestions to improve these codes. The model I created was a 1D radial flow to a well. Also, any idea on how to use other optimization algorithms in conjunction to the mph to iterative update the parameters to minimize the likelihood function will be greatly appreciated. The codes are below.

Import necessary libraries

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import minimize
import mph

=======================

1. COMSOL INTEGRATION

=======================

Launch MPH using the license

mph.option('classkit', True)
client = mph.start()

Load the COMSOL model

model = client.load(r"D:\COMSOL\Parameter_estimation1.mph")

Solve the model

model.solve()

=======================

2. LOAD OBSERVED DATA

=======================

Observed field data (head values at different distances)

observed_data = pd.DataFrame({
"Distance (m)": [0.10, 66.76, 133.42, 200.08, 266.74, 333.40, 400.06, 466.72, 533.38, 600.04,
666.70, 733.36, 800.02, 866.68, 933.34, 1000.00],
"Head (m)": [-0.009053, -0.004306, -0.003206, -0.002561, -0.002103, -0.001748, -0.001458,
-0.001213, -0.001000, -0.000813, -0.000645, -0.000494, -0.000355, -0.000228,
-0.000110, 0.000000]
})

Extract distances for simulation

locations = observed_data["Distance (m)"].values
observed_heads = observed_data["Head (m)"].values

=======================

3. DEFINE PRIOR DISTRIBUTIONS

=======================

Prior assumptions based on given conditions

prior_k = 1e-4 # Hydraulic conductivity (m/s)
prior_well_flux = -0.001 # Well flux (m³/s)
prior_aquifer_thickness = 10 # Aquifer thickness (m)
prior_radius_influence = 1000 # Radius of influence (m)
prior_well_radius = 0.1 # Initial well radius (m)

=======================

4. SIMULATION FUNCTION

=======================

def simulate_model(k):
"""
Run COMSOL simulation for a given hydraulic conductivity 'k'
and return the simulated head values at specified locations.
"""
# Update Hydraulic Conductivity in COMSOL model
model.parameter('K', f'{k} [m/s]')

# Run simulation
model.solve()

# Extract solution dataset (Ensure correct node name in COMSOL)
solution = model.evaluate('Data', dataset='Study 1//Solution 1')

# The comsol results always have some headings - need to remove all them or make the algorithm automatically detects only the data points.
# Ensure the output is properly structured
if len(solution) != len(locations):
    raise ValueError("Mismatch between simulation output and observation points.")

return np.array(solution)

=======================

5. LIKELIHOOD FUNCTION

=======================

def likelihood(k, observed_heads, sigma=0.01):
"""
Compute the likelihood of a given hydraulic conductivity 'k'
based on observed vs. simulated head values.
"""
try:
simulated_heads = simulate_model(k)
except Exception as e:
print(f"Error in COMSOL simulation: {e}")
return -np.inf # Assign very low likelihood to avoid errors

# Compute residuals
residuals = observed_heads - simulated_heads

# Compute log-likelihood assuming normal errors
log_likelihood = np.sum(norm.logpdf(residuals, loc=0, scale=sigma))

return log_likelihood

=======================

6. MCMC PARAMETER ESTIMATION

=======================

def mcmc_parameter_estimation(initial_k, observed_heads, n_iterations=5000, step_size=1e-5):
"""
Bayesian MCMC using Metropolis-Hastings algorithm for parameter estimation.
"""
k_samples = []
current_k = initial_k
current_likelihood = likelihood(current_k, observed_heads)

for i in range(n_iterations):
    # Propose new value for k using Gaussian random walk
    proposed_k = np.abs(np.random.normal(current_k, step_size))

    # Compute likelihood for proposed k
    proposed_likelihood = likelihood(proposed_k, observed_heads)

    # Compute acceptance probability
    acceptance_ratio = np.exp(proposed_likelihood - current_likelihood)

    # Accept/reject step
    if np.random.rand() < acceptance_ratio:
        current_k = proposed_k
        current_likelihood = proposed_likelihood

    k_samples.append(current_k)

    # Print progress every 500 iterations
    if i % 500 == 0:
        print(f"Iteration {i}: Current K = {current_k:.6e}, Likelihood = {current_likelihood:.6e}")

return np.array(k_samples)

=======================

7. RUN MCMC

=======================

initial_k = 1e-5 # Initial guess for hydraulic conductivity
k_samples = mcmc_parameter_estimation(initial_k, observed_heads)

=======================

8. VISUALIZE RESULTS

=======================

Plot observed vs. simulated fit using mean K

plt.figure(figsize=(8, 6))
plt.scatter(locations, observed_heads, label="Observed Data", color='red')
plt.plot(locations, simulate_model(np.mean(k_samples)), label="Simulated Fit (Mean K)", color='blue')
plt.xlabel("Distance (m)")
plt.ylabel("Head (m)")
plt.legend()
plt.title("Hydraulic Conductivity Estimation")
plt.show()

Plot posterior distribution of K

plt.figure(figsize=(8, 6))
plt.hist(k_samples, bins=30, density=True)
plt.xlabel("Hydraulic Conductivity (m/s)")
plt.ylabel("Density")
plt.title("MCMC Posterior Distribution of Hydraulic Conductivity")
plt.show()