Python实现油气体积计量蒙特卡洛分析的技术求助

Hey Josiah, let's walk through fixing your code and building a straightforward Monte Carlo simulation tailored to your reservoir volume calculation. Since you're coming from VBA, we'll keep things clear and avoid overly complex Python tricks.

First: What's Wrong with Your Current Code? #

Let's break down the key issues:

• Import Syntax Error: You can't combine two imports in one line like from scipy.stats import * import numpy as np—they need separate lines.
• No Random Samples: norm(200,50) creates a distribution object, but doesn't generate the 10000 random values you need for simulation. You need the .rvs() method for that.
• Useless Loop: Your for loop defines distributions but doesn't run any calculations, and the return statement breaks the loop immediately.
• Invalid Calculations: You're trying to do math with distribution objects instead of numerical arrays, which will throw errors.

Working Beginner-Friendly Code #

Here's a corrected version with comments explaining every step:

# Import required libraries: numpy for arrays/statistics, norm for normal distributions
import numpy as np
from scipy.stats import norm

# Set number of Monte Carlo simulations
n_simulations = 10000

# Generate 10000 random samples for each variable (normal distribution)
# loc = mean, scale = standard deviation
Area = norm(loc=200, scale=50).rvs(n_simulations)
Thickness = norm(loc=100, scale=25).rvs(n_simulations)
# Clip NTG to 0-1 since it can't be negative or exceed 1
NTG = np.clip(norm(loc=0.85, scale=0.2).rvs(n_simulations), 0, 1)
# Clip POR to 0-1 for physical plausibility
POR = np.clip(norm(loc=0.32, scale=0.02).rvs(n_simulations), 0, 1)
GS = np.clip(norm(loc=0.80, scale=0.2).rvs(n_simulations), 0, 1)
BG = norm(loc=0.0024, scale=0.0001).rvs(n_simulations)
# Feather is a constant, so create an array of 1s with 10000 elements
Feather = np.ones(n_simulations)

# Calculate volume for ALL simulations at once (vectorized operation, no loops needed!)
# This is way faster than looping in VBA
volume = (((Area * Thickness * NTG * POR * GS) / BG) * 43560 * Feather) / 1000000000

# Calculate P90, P50, P10 percentiles
# Note: P90 = 10th percentile (90% of results are >= this value)
# P50 = median, P10 = 90th percentile (10% of results are >= this value)
p90 = np.percentile(volume, 10)
p50 = np.percentile(volume, 50)
p10 = np.percentile(volume, 90)

# Print results with clean formatting
print(f"Monte Carlo Simulation Results ({n_simulations} runs):")
print(f"P90 Volume: {p90:.2f} billion cubic feet")
print(f"P50 Volume: {p50:.2f} billion cubic feet")
print(f"P10 Volume: {p10:.2f} billion cubic feet")

Key Concepts Explained for VBA Users #

• Vectorized Operations: Instead of writing a For loop to calculate each simulation one by one, numpy lets you do math on entire arrays at once. This is way faster and cleaner than VBA loops.
• Random Samples: .rvs(n_simulations) generates n_simulations random numbers from the normal distribution you defined.
• Clipping Values: Variables like NTG (net-to-gross) can't be negative or greater than 1. np.clip() ensures your samples stay within physically realistic bounds.
• Percentiles: np.percentile() calculates the value where X% of the data falls below it. For reservoir engineering:
  • P90: Conservative estimate (90% chance the actual volume is at least this value)
  • P50: Most likely estimate (median)
  • P10: Optimistic estimate (10% chance the actual volume is at least this value)

Extra Tip for Validation #

If you want to check the distribution of your results, you can add a quick histogram (great for debugging):

import matplotlib.pyplot as plt

plt.hist(volume, bins=30, edgecolor='black')
plt.title('Distribution of Simulated Reservoir Volumes')
plt.xlabel('Volume (billion cubic feet)')
plt.ylabel('Number of Simulations')
plt.axvline(p90, color='red', linestyle='--', label='P90')
plt.axvline(p50, color='green', linestyle='--', label='P50')
plt.axvline(p10, color='blue', linestyle='--', label='P10')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Josiah Hulsey