Portfolio DCA

Level 3

Chapter 9: Decline Curve Analysis

descriptionProblem

Generate synthetic production data for 20 wells with varied parameters. Run the auto_dca() function on all of them. Create a summary report that includes: (a) a bar chart ranking wells by remaining reserves, (b) a histogram of $b$ -factors across the portfolio, and (c) the total field EUR and cumulative production.

---

You run reserves for the OML 58 development. There are 20 producing wells (OD-001 … OD-020), and reserves review is next week. Doing decline analysis by hand on each well does not scale, so you build a portfolio engine: one function that fits any well, and one that rolls all 20 up into a ranked reserves report.

## The portfolio

Twenty deterministic synthetic wells are already defined for you in WELLS. Each is 36 months of noise-free hyperbolic production generated from known Arps parameters (well i, i = 0 … 19):

qi_bopd = 1000 + 50*i
di_per_month = 0.08 + 0.004*i
b_factor = min(0.4 + 0.02*i, 0.95)

Each entry is a dict {"name", "t_months", "q_bopd"} where t_months is 0, 1, …, 35 and q_bopd is the hyperbolic rate from arps_hyperbolic.

## What to build

auto_dca(t_months, q_bopd, q_abandon_bopd=20.0) → dict. Fit the

hyperbolic model to one well with scipy.optimize.curve_fit, then compute its EUR with eur_hyperbolic. Return:

key	meaning
`qi_bopd`	fitted initial rate
`di_per_month`	fitted initial decline
`b_factor`	fitted Arps b-exponent
`eur_bbl`	`eur_hyperbolic(qi, Di, b, q_abandon) * 30.44`

Rates qi_bopd are in bbl/day and t/Di are per month, so the raw eur_hyperbolic result is in bopd·months. Multiply by DAYS_PER_MONTH = 30.44 to report eur_bbl as a true volume in barrels (this matches the book's day/month convention).

Use bounds that keep the fit physical: Di > 0 and 0 < b < 1 (so the hyperbolic EUR stays finite). A good start point is p0 = [q_bopd[0], 0.1, 0.5].

portfolio(wells) → dict. Run auto_dca on every well and roll up:

key	meaning
`summary`	list of dicts, one per well, each the `auto_dca` result plus the well's `name` (same order as `wells`)
`total_eur_bbl`	sum of every well's `eur_bbl`
`ranking`	the same per-well dicts sorted by `eur_bbl` descending

Run it on WELLS and store the result in result.

Two plots (use plt.subplots(1, 2, ...)): a bar chart of

eur_bbl per well, and a histogram of the fitted b_factor across the portfolio. Label both axes on both panels.

Because the data is noise-free, curve_fit recovers the true parameters almost exactly: your fitted b_factor for well 0 should land near 0.40, and the whole field EUR should total roughly 16 million barrels (after the * 30.44 day/month conversion).

lightbulbHints (0/3)

Stuck? Reveal hints one at a time — they progress from nudge to near-solution.

codeYour solution

main.py

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# ── Embedded Arps decline functions (do not edit) ────────────────────────
def arps_hyperbolic(t, qi, Di, b):
    return qi / (1.0 + b * Di * t) ** (1.0 / b)

def eur_hyperbolic(qi, Di, b, q_abandon):
    if qi <= q_abandon or b >= 1:
        return float("nan")
    return (qi ** b / (Di * (1.0 - b))) * (qi ** (1.0 - b) - q_abandon ** (1.0 - b))

# ── OML 58 portfolio: 20 deterministic noise-free hyperbolic wells ───────
# Well i (i = 0..19): qi = 1000+50i, Di = 0.08+0.004i, b = min(0.4+0.02i, 0.95).
def _build_wells():
    t_months = np.arange(36, dtype=float)
    wells = []
    for i in range(20):
        qi_bopd = 1000.0 + 50.0 * i
        di_per_month = 0.08 + 0.004 * i
        b_factor = min(0.4 + 0.02 * i, 0.95)
        q_bopd = arps_hyperbolic(t_months, qi_bopd, di_per_month, b_factor)
        wells.append({
            "name": f"OD-{i + 1:03d}",
            "t_months": t_months,
            "q_bopd": q_bopd,
        })
    return wells

WELLS = _build_wells()
Q_ABANDON_BOPD = 20.0
# q is in bbl/DAY but t and Di are per-MONTH, so eur_hyperbolic returns
# bopd*months - multiply by days/month to report a true barrel volume.
DAYS_PER_MONTH = 30.44

def auto_dca(t_months, q_bopd, q_abandon_bopd=Q_ABANDON_BOPD):
    """Fit hyperbolic Arps to one well, return {qi_bopd, di_per_month, b_factor, eur_bbl}."""
    # TODO: curve_fit arps_hyperbolic to (t_months, q_bopd)
    #       p0 = [q_bopd[0], 0.1, 0.5]; bounds keep Di>0 and 0<b<1.
    # TODO: unpack fitted qi, Di, b
    # TODO: eur_bbl = eur_hyperbolic(qi, Di, b, q_abandon_bopd) * DAYS_PER_MONTH
    #       (eur_hyperbolic returns bopd*months; * DAYS_PER_MONTH -> barrels)
    # TODO: return the four-key dict
    pass

def portfolio(wells):
    """Run auto_dca on every well; return {summary, total_eur_bbl, ranking}."""
    # TODO: summary = [auto_dca(...) + {"name": ...} for each well] (input order)
    # TODO: total_eur_bbl = sum of every well's eur_bbl
    # TODO: ranking = summary sorted by eur_bbl descending
    pass

result = None  # TODO: portfolio(WELLS)

# TODO: fig, axes = plt.subplots(1, 2, figsize=(14, 5))
#       axes[0] -> bar chart of eur_bbl per well
#       axes[1] -> histogram of fitted b_factor
#       label both axes on both panels; plt.show()

print("wells:", len(WELLS))
print("portfolio total EUR (bbl):",
      None if result is None else result["total_eur_bbl"])
print("top well:",
      None if result is None else result["ranking"][0]["name"])

visibilityReveal reference solutionexpand_more

Try solving it yourself first — the hints walk you through it. The solution below is one valid approach; yours may differ and still be correct.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


# ── Embedded Arps decline functions (do not edit) ────────────────────────
def arps_hyperbolic(t, qi, Di, b):
    return qi / (1.0 + b * Di * t) ** (1.0 / b)


def eur_hyperbolic(qi, Di, b, q_abandon):
    if qi <= q_abandon or b >= 1:
        return float("nan")
    return (qi ** b / (Di * (1.0 - b))) * (qi ** (1.0 - b) - q_abandon ** (1.0 - b))


# ── OML 58 portfolio: 20 deterministic noise-free hyperbolic wells ───────
# Well i (i = 0..19): qi = 1000+50i, Di = 0.08+0.004i, b = min(0.4+0.02i, 0.95).
def _build_wells():
    t_months = np.arange(36, dtype=float)
    wells = []
    for i in range(20):
        qi_bopd = 1000.0 + 50.0 * i
        di_per_month = 0.08 + 0.004 * i
        b_factor = min(0.4 + 0.02 * i, 0.95)
        q_bopd = arps_hyperbolic(t_months, qi_bopd, di_per_month, b_factor)
        wells.append({
            "name": f"OD-{i + 1:03d}",
            "t_months": t_months,
            "q_bopd": q_bopd,
        })
    return wells


WELLS = _build_wells()
Q_ABANDON_BOPD = 20.0
# q is in bbl/DAY but t and Di are per-MONTH, so eur_hyperbolic returns
# bopd*months - multiply by days/month to report a true barrel volume.
DAYS_PER_MONTH = 30.44


def auto_dca(t_months, q_bopd, q_abandon_bopd=Q_ABANDON_BOPD):
    """Fit hyperbolic Arps to one well, return {qi_bopd, di_per_month, b_factor, eur_bbl}."""
    t = np.asarray(t_months, dtype=float)
    q = np.asarray(q_bopd, dtype=float)
    p0 = [float(q[0]), 0.1, 0.5]
    bounds = ([1e-6, 1e-6, 1e-6], [np.inf, 5.0, 0.999])
    popt, _ = curve_fit(arps_hyperbolic, t, q, p0=p0, bounds=bounds, maxfev=100000)
    qi_bopd = float(popt[0])
    di_per_month = float(popt[1])
    b_factor = float(popt[2])
    # eur_hyperbolic returns bopd*months; * DAYS_PER_MONTH -> barrels.
    eur_bbl = float(eur_hyperbolic(qi_bopd, di_per_month, b_factor, q_abandon_bopd)
                    * DAYS_PER_MONTH)
    return {
        "qi_bopd": qi_bopd,
        "di_per_month": di_per_month,
        "b_factor": b_factor,
        "eur_bbl": eur_bbl,
    }


def portfolio(wells):
    """Run auto_dca on every well; return {summary, total_eur_bbl, ranking}."""
    summary = []
    for w in wells:
        fit = auto_dca(w["t_months"], w["q_bopd"])
        fit = dict(fit)
        fit["name"] = w["name"]
        summary.append(fit)
    total_eur_bbl = float(sum(row["eur_bbl"] for row in summary))
    ranking = sorted(summary, key=lambda row: row["eur_bbl"], reverse=True)
    return {
        "summary": summary,
        "total_eur_bbl": total_eur_bbl,
        "ranking": ranking,
    }


result = portfolio(WELLS)

# ── Portfolio plots: EUR bar chart + b-factor histogram ──────────────────
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

names = [row["name"] for row in result["summary"]]
eurs = [row["eur_bbl"] for row in result["summary"]]
axes[0].bar(range(len(names)), eurs)
axes[0].set_xlabel("Well")
axes[0].set_ylabel("EUR (bbl)")
axes[0].set_title("EUR by well - OML 58 portfolio")

b_factors = [row["b_factor"] for row in result["summary"]]
axes[1].hist(b_factors, bins=8)
axes[1].set_xlabel("Arps b-factor")
axes[1].set_ylabel("Well count")
axes[1].set_title("b-factor distribution")

fig.suptitle("OML 58 portfolio decline analysis", fontweight="bold")
plt.tight_layout()
plt.show()

print("wells:", len(WELLS))
print("portfolio total EUR (bbl):", result["total_eur_bbl"])
print("top well:", result["ranking"][0]["name"])

lockCopying code is a Full Access feature.

arrow_back

9.7 Detecting a Decline-Regime Change

9.9 DCA vs Volumetric Reserves

arrow_forward