Ex.2 演習ノート ―Volveフィールドデータで代理モデルを構築する―¶

実践 貯留層工学 データ駆動型シミュレータ編 第2回 演習(Ex.2)

本ノートは データ駆動型シミュレータ編_第2回_演習.pptx(パートF)の構成 (データ再読込→入力・出力設計→RSM/回帰→ランダムフォレスト→(任意)NN→精度検証→時系列分割の検証)に沿って、 第1回演習(データ駆動型シミュレータ編_第1回_演習notebook.ipynb)で読み込んだVolveフィールドデータのうち 油田全体の月次生産量・圧力データを用いて、実際に代理モデルを構築・検証した記録である。

使用データ: Volve_Production_Data_Processed.csv(月次112行、2007年9月〜2016年12月)

Step1: 環境構築確認・第1回からの接続¶

In [3]:
import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import sklearn
import warnings

from sklearn.model_selection import train_test_split, KFold, cross_val_score
from sklearn.preprocessing import StandardScaler, PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.neural_network import MLPRegressor
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
from sklearn.pipeline import make_pipeline

plt.rcParams["font.family"] = "Yu Gothic"
plt.rcParams["axes.unicode_minus"] = False

print("pandas", pd.__version__, "/ numpy", np.__version__,
      "/ matplotlib", matplotlib.__version__, "/ scikit-learn", sklearn.__version__)

DATA_DIR = r"C:\kawata\testAI\調査部門\参考資料\Volve"

pandas 3.0.3 / numpy 2.5.0 / matplotlib 3.11.1 / scikit-learn 1.9.0

Step2: データ再読込・特徴量設計¶

第1回Ex.1では坑井別の坑底圧力データ(Volve_BHFP_data.csv)・坑井別特徴量テーブル(7坑井)を扱った。 7坑井分の集計テーブルでは代理モデルの学習データとして行数(サンプル数)が少なすぎるため、 第2回Ex.2では油田全体の月次生産量・圧力データ(Volve_Production_Data_Processed.csv、112行)を主教材に切り替える。

前処理の注意(第1回と同じ落とし穴): 本ファイルも日付が DD/MM/YYYY 表記のため、 第1回で気づいた dayfirst=True を引き続き指定する必要がある。

In [5]:
df = pd.read_csv(f"{DATA_DIR}/Volve_Production_Data_Processed.csv")
df["Date"] = pd.to_datetime(df["Date"], dayfirst=True)  # 第1回の気づき: dayfirst=True が必須
df["elapsed_days"] = (df["Date"] - df["Date"].min()).dt.days
print(df.shape)
print(df["Date"].min(), "〜", df["Date"].max())
df.head()

(112, 9)
2007-09-01 00:00:00 〜 2016-12-01 00:00:00

Step3: 入力・出力変数の設計と多重共線性の確認¶

  • 出力変数(y): p (psia) ―貯留層圧力
  • 入力変数(X): elapsed_days(経過日数)、Np (STB)・Gp (SCF)・Wp (STB)(累積生産量)、Gi (SCF)・Wi (STB)(累積圧入量)
In [7]:
feature_cols = ["elapsed_days", "Np (STB)", "Gp (SCF)", "Wp (STB)", "Gi (SCF)", "Wi (STB)"]
X = df[feature_cols].values
y = df["p (psia)"].values

print("Gi(累積ガス圧入量)が非ゼロの行数:", (df["Gi (SCF)"] != 0).sum(), "/", len(df))
corr = df[["Np (STB)", "Gp (SCF)", "Wp (STB)", "Wi (STB)"]].corr()
print(corr.round(4))

Gi(累積ガス圧入量)が非ゼロの行数: 0 / 112
          Np (STB)  Gp (SCF)  Wp (STB)  Wi (STB)
Np (STB)    1.0000    0.9998    0.8814    0.9664
Gp (SCF)    0.9998    1.0000    0.8902    0.9711
Wp (STB)    0.8814    0.8902    1.0000    0.9731
Wi (STB)    0.9664    0.9711    0.9731    1.0000

気づき(多重共線性): Np (STB) と Gp (SCF) の相関係数は0.9998とほぼ完全な共線関係にある (Volveは溶存ガスを主体とするため、油量が決まるとガス量もほぼ一意に決まる)。 Gi (SCF)(累積ガス圧入量)は全112行で常に0 ―Volveは水攻(water injection)主体でガス圧入を行っていないためで、 この変数は代理モデルの入力としては意味を持たない(Step5の特徴量重要度でも重要度0として現れる)。 パートB6で扱った「入力変数の中心化・標準化により多重共線性を軽減する」という指摘を、実データで確認できた。

Step4: RSM(多項式回帰)・線形回帰モデルの構築(ホールドアウト検証)¶

In [10]:
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
print("訓練データ:", len(X_train), "件 / テストデータ:", len(X_test), "件")

def evaluate(name, model, X_tr, y_tr, X_te, y_te):
    model.fit(X_tr, y_tr)
    pred = model.predict(X_te)
    rmse = mean_squared_error(y_te, pred) ** 0.5
    mae = mean_absolute_error(y_te, pred)
    r2 = r2_score(y_te, pred)
    print(f"{name:10s}  RMSE={rmse:7.2f} psia  MAE={mae:7.2f} psia  R2={r2:6.3f}")
    return dict(name=name, model=model, rmse=rmse, mae=mae, r2=r2, pred=pred)

results = {}
results["linear"] = evaluate("線形回帰", make_pipeline(StandardScaler(), LinearRegression()),
                               X_train, y_train, X_test, y_test)
results["poly2"] = evaluate("多項式(2次)", make_pipeline(StandardScaler(), PolynomialFeatures(degree=2, include_bias=False), LinearRegression()),
                              X_train, y_train, X_test, y_test)

訓練データ: 89 件 / テストデータ: 23 件
線形回帰        RMSE=  89.74 psia  MAE=  51.21 psia  R2= 0.831
多項式(2次)     RMSE=  75.51 psia  MAE=  48.77 psia  R2= 0.881

Step5: ランダムフォレストモデルの構築(ホールドアウト検証)¶

In [12]:
rf = RandomForestRegressor(n_estimators=300, random_state=42, oob_score=True)
results["rf"] = evaluate("RF", rf, X_train, y_train, X_test, y_test)
print("OOBスコア:", round(rf.oob_score_, 3))

importances = pd.Series(rf.feature_importances_, index=feature_cols).sort_values(ascending=False)
print("\n特徴量重要度:")
print(importances.round(3))

RF          RMSE=  49.32 psia  MAE=  31.11 psia  R2= 0.949
OOBスコア: 0.91

特徴量重要度:
elapsed_days    0.231
Np (STB)        0.216
Gp (SCF)        0.191
Wi (STB)        0.182
Wp (STB)        0.180
Gi (SCF)        0.000
dtype: float64

気づき: ランダムフォレストが最も高精度(R2=0.949、RMSE≈49psia)。 特徴量重要度は経過日数・Np・Gp・Wi・Wpが拮抗する一方、Gi (SCF)は重要度0.000 ―Step3で確認した 「全期間ゼロの変数」という実データの構造が、そのままモデルの重要度にも反映されている。

Step6: ニューラルネットワークモデルの構築(任意課題・ホールドアウト検証)¶

In [15]:
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    mlp = make_pipeline(StandardScaler(), MLPRegressor(hidden_layer_sizes=(32, 16), max_iter=5000, random_state=42))
    results["mlp"] = evaluate("MLP(NN)", mlp, X_train, y_train, X_test, y_test)

MLP(NN)     RMSE= 137.20 psia  MAE= 109.06 psia  R2= 0.606

気づき: ニューラルネットワーク(MLP)はRMSE≈137psia・R2=0.606と、RSM(R2=0.881)・ ランダムフォレスト(R2=0.949)を下回った。訓練データが89件と少なく、パートC10で扱った 「ニューラルネットワークは大規模データを要し、中小規模データではランダムフォレストの方が頑健」という 一般論を、実データで裏付ける結果になった。

Step7: 5-fold交差検証(ランダムフォレスト)¶

In [18]:
kf = KFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(RandomForestRegressor(n_estimators=300, random_state=42), X, y, cv=kf, scoring="r2")
print("5-fold R2:", np.round(scores, 3))
print("平均R2:", round(scores.mean(), 3), " / 標準偏差:", round(scores.std(), 3))

5-fold R2: [0.95  0.901 0.768 0.897 0.972]
平均R2: 0.897  / 標準偏差: 0.071

気づき: ホールドアウト法単独ではR2=0.949だったが、5-fold交差検証ではfoldごとにR2が 0.768〜0.972まで変動し、平均は0.897だった。単一の分割結果だけに頼らず、複数分割での確認が 精度評価の頑健性を高めることを実データで確認できた(パートE2の内容の実例)。

Step8: 時系列分割での検証(重要な気づき)¶

ここまではランダムな80/20分割で検証したが、Volveのデータは時系列データである。 実務での予測(将来の圧力を予測する用途)によりふさわしい検証として、時系列に沿って 前半80%を訓練・残り20%(2015年2月以降)をテストとする分割も試す。

In [21]:
n = len(df)
cut = int(n * 0.8)
X_train_t, X_test_t = X[:cut], X[cut:]
y_train_t, y_test_t = y[:cut], y[cut:]
print("訓練期間:", df["Date"].iloc[0].date(), "〜", df["Date"].iloc[cut-1].date())
print("テスト期間:", df["Date"].iloc[cut].date(), "〜", df["Date"].iloc[-1].date())
print("訓練データの圧力範囲: {:.2f} 〜 {:.2f} psia".format(y_train_t.min(), y_train_t.max()))
print("テストデータの圧力範囲: {:.2f} 〜 {:.2f} psia".format(y_test_t.min(), y_test_t.max()))

rf_t = RandomForestRegressor(n_estimators=300, random_state=42)
poly_t = make_pipeline(StandardScaler(), PolynomialFeatures(degree=2, include_bias=False), LinearRegression())

time_results = {}
for name, model in [("RF", rf_t), ("多項式(2次)", poly_t)]:
    model.fit(X_train_t, y_train_t)
    pred = model.predict(X_test_t)
    rmse = mean_squared_error(y_test_t, pred) ** 0.5
    r2 = r2_score(y_test_t, pred)
    time_results[name] = dict(pred=pred, rmse=rmse, r2=r2)
    print(f"[時系列分割] {name:10s}  RMSE={rmse:7.2f} psia  R2={r2:7.3f}")

訓練期間: 2007-09-01 〜 2015-01-01
テスト期間: 2015-02-01 〜 2016-12-01
訓練データの圧力範囲: 4085.44 〜 5114.90 psia
テストデータの圧力範囲: 4815.37 〜 5273.87 psia
[時系列分割] RF          RMSE= 134.26 psia  R2= -0.419
[時系列分割] 多項式(2次)     RMSE= 408.79 psia  R2=-12.153

気づき(重要): ランダム分割ではRF R2=0.949と良好だったが、時系列分割ではRF R2=-0.42、 多項式回帰はR2=-12.15まで悪化した。原因は、テスト期間(2015年2月〜2016年12月)の圧力が 最大5273.87psiaに達しているのに対し、訓練期間(2007年9月〜2015年1月)の圧力は最大5114.90psiaまでしか 経験しておらず、モデルが**学習データ範囲外(外挿領域)**での予測を強いられているためである。 パートE5・B10で「学習データの範囲外への適用には注意」と概念として説明していた内容を、 Volveの実データで定量的に確認できた。ランダム分割だけで精度を評価すると、この種の外挿リスクを 見落とす恐れがある、という本演習最大の教訓である。

Step9: 結果の可視化¶

In [24]:
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
colors = {"linear": "#4A3AA7", "poly2": "#2A78D6", "rf": "#1BAF7A", "mlp": "#EB6834"}
labels = {"linear": "線形回帰", "poly2": "多項式(2次)", "rf": "ランダムフォレスト", "mlp": "NN(MLP)"}
for key, res in results.items():
    ax.scatter(y_test, res["pred"], s=28, alpha=0.75, color=colors[key], label=f"{labels[key]} (R2={res['r2']:.2f})")
lims = [min(y_test.min(), y_test.min()) - 50, y_test.max() + 50]
ax.plot(lims, lims, color="#89877E", linestyle="--", linewidth=1)
ax.set_xlabel("実測値 p (psia)"); ax.set_ylabel("予測値 p (psia)")
ax.set_title("ランダム分割: 予測値 vs 実測値")
ax.legend(fontsize=8); ax.grid(alpha=0.3)

ax2 = axes[1]
names = list(results.keys())
rmses = [results[k]["rmse"] for k in names]
ax2.bar([labels[k] for k in names], rmses, color=[colors[k] for k in names])
ax2.set_ylabel("RMSE (psia)")
ax2.set_title("手法別のRMSE比較(ランダム分割)")
for i, v in enumerate(rmses):
    ax2.text(i, v + 2, f"{v:.1f}", ha="center", fontsize=9)
plt.setp(ax2.get_xticklabels(), rotation=15, ha="right")

fig.tight_layout()
plt.show()
No description has been provided for this image
In [25]:
fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(df["Date"], df["p (psia)"], color="#89877E", linewidth=1.2, label="実測値 p (psia)")
ax.plot(df["Date"].iloc[cut:], time_results["RF"]["pred"], color="#1BAF7A", linewidth=1.6,
        linestyle="--", label="RF予測(時系列分割テスト区間)")
ax.axvline(df["Date"].iloc[cut], color="#EB6834", linestyle=":", linewidth=1.5)
ax.text(df["Date"].iloc[cut], ax.get_ylim()[1], " ← 訓練/テスト境界", color="#EB6834", fontsize=9, va="top")
ax.set_xlabel("日付"); ax.set_ylabel("貯留層圧力 p (psia)")
ax.set_title("時系列分割: テスト期間で実測値から外挿し予測が乖離する様子")
ax.legend(fontsize=9); ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()
No description has been provided for this image

まとめ・第3回への申し送り¶

  • Volveの月次生産量・圧力データ(112行)を用いて、RSM(線形・2次多項式回帰)・ランダムフォレスト・ ニューラルネットワーク(任意課題)の3手法で貯留層圧力pの代理モデルを構築し、ホールドアウト法・ 5-fold交差検証で精度を検証する一連のワークフローを、pptx資料(F1〜F10)の流れに沿って実行できた。
  • 実データ検証で分かった3点。第2回の教材化・受講者への注意喚起に反映すべき:
    1. 坑井別特徴量テーブル(7坑井)はサンプル数が少なすぎるため、月次112行の油田全体データを 代理モデル構築の主教材とした。
    2. Np (STB)とGp (SCF)は相関係数0.9998とほぼ共線関係にあり、Gi (SCF)は全期間ゼロ (Volveは水攻主体でガス圧入を行っていない)―多重共線性・無情報変数の実例として教材に使える。
    3. **ランダム分割ではRF R2=0.949と良好だが、時系列分割(2015年2月以降をテスト)ではRF R2=-0.42、 多項式回帰はR2=-12.15まで悪化する。**原因はテスト期間の圧力が訓練期間の範囲を超える外挿。 検証方法(ランダム分割 vs 時系列分割)次第で結論が一変するという、本演習最大の教訓が得られた。
  • この「代理モデルの信頼できる範囲」「検証設計次第で結果が変わる」という気づきは、第3回 (データ駆動型ヒストリーマッチング)での不確実性評価・予測範囲の議論にそのまま接続する。