Production Optimization

schedule15 min readfitness_center10 exercises

Why This Chapter Exists

A petroleum engineer responsible for a producing field faces a daily allocation problem. Multiple wells share limited processing capacity --- gas handling, water treatment, pipeline throughput, power supply. Each well has different characteristics: different reservoir pressures, different water cuts, different gas-oil ratios, different decline rates. The question is not simply "how much can each well produce?" but "how should production be distributed across all wells to maximize total output within the constraints the facilities impose?"

This is an optimization problem, and it is one of the most practically valuable skills a production engineer can have. A field producing 20,000 barrels per day under a naive allocation might produce 23,000 under an optimized one --- with the same wells, the same facilities, and the same reservoir. The difference is entirely in how the rates are distributed.

This chapter teaches you to formulate production problems as mathematical optimization problems, solve them in Python, and interpret the results in engineering terms. You will also build automated surveillance tools that flag when a well is underperforming relative to its potential.

infoWhat You Will Learn

Formulate production allocation as a constrained optimization problem
Use scipy.optimize for linear and nonlinear optimization
Allocate rates across wells to maximize total oil production
Optimize gas lift injection rates across a field
Build production surveillance dashboards
Detect underperforming wells using statistical methods

The Allocation Problem

Consider a field with five producing wells feeding into a single processing facility. The facility can handle at most 3,000 barrels of total liquid per day, 5 MMscf of gas per day, and 1,200 barrels of water per day. Each well has a different IPR, a different gas-oil ratio, and a different water cut. Some wells are more oil-rich but gas-heavy. Others produce mostly water with modest oil.

The naive approach is to let each well produce at its natural rate. But the natural rates may exceed the facility limits, requiring some wells to be choked back. The question becomes: which wells do you choke, and by how much?

The optimal answer is not obvious. Choking the highest-GOR well saves gas capacity but may sacrifice more oil than choking a lower-rate well with moderate GOR. The right allocation depends on the interaction between all the constraints simultaneously.

main.py

import numpy as np
import pandas as pd
from scipy.optimize import linprog, minimize

# Field data: 5 wells with different characteristics
wells = pd.DataFrame({
    "Well": ["A-01", "A-02", "A-03", "B-01", "B-02"],
    "Max Oil Rate (STB/d)": [1200, 800, 600, 1500, 900],
    "GOR (scf/STB)": [450, 1200, 300, 800, 600],
    "Water Cut (frac)": [0.15, 0.05, 0.45, 0.10, 0.30],
    "Operating Cost ($/STB)": [12, 18, 8, 15, 10],
})

# Facility constraints
max_liquid_bpd = 3000
max_gas_mmscfd = 5.0    # 5,000 Mscf/d
max_water_bpd = 1200

print("FIELD DATA")
print(wells.to_string(index=False))
print()
print("FACILITY CONSTRAINTS")
print(f"  Max total liquid:  {max_liquid_bpd:,} bbl/day")
print(f"  Max gas handling:  {max_gas_mmscfd} MMscf/day")
print(f"  Max water:         {max_water_bpd:,} bbl/day")

Linear Optimization — Maximizing Oil Production

When the relationship between the decision variables (well rates) and the objective (total oil) is linear, and the constraints are linear inequalities, the problem can be solved using linear programming. This is the simplest and most common formulation for production allocation.

The objective is to maximize total oil production:

\text{Maximize: } \sum_{i=1}^{n} q_{o,i}

subject to:

\sum_{i=1}^{n} q_{o,i} \times \frac{1}{1 - WC_i} \leq Q_{liquid,max} \quad \text{(total liquid constraint)}

\sum_{i=1}^{n} q_{o,i} \times GOR_i \leq Q_{gas,max} \quad \text{(gas handling constraint)}

\sum_{i=1}^{n} q_{o,i} \times \frac{WC_i}{1 - WC_i} \leq Q_{water,max} \quad \text{(water disposal constraint)}

0 \leq q_{o,i} \leq q_{o,i,max} \quad \text{(individual well limits)}

main.py

import numpy as np
import pandas as pd
from scipy.optimize import linprog

# Extract well data
max_oil = wells["Max Oil Rate (STB/d)"].values
gor = wells["GOR (scf/STB)"].values
wc = wells["Water Cut (frac)"].values
n_wells = len(max_oil)

# Objective: maximize total oil = minimize negative total oil
c = -np.ones(n_wells)  # linprog minimizes, so negate

# Inequality constraints: A_ub @ x <= b_ub
# 1. Total liquid: qo / (1 - WC) for each well
liquid_coeffs = 1.0 / (1.0 - wc)

# 2. Gas handling: qo * GOR (convert scf to MMscf)
gas_coeffs = gor / 1e6

# 3. Water: qo * WC / (1 - WC)
water_coeffs = wc / (1.0 - wc)

A_ub = np.array([liquid_coeffs, gas_coeffs, water_coeffs])
b_ub = np.array([max_liquid_bpd, max_gas_mmscfd, max_water_bpd])

# Bounds: 0 <= qo <= max_oil for each well
bounds = [(0, m) for m in max_oil]

# Solve
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if result.success:
    optimal_rates = result.x
    total_oil = sum(optimal_rates)

    # Build results table
    results = pd.DataFrame({
        "Well": wells["Well"],
        "Max Rate": max_oil,
        "Optimal Rate": np.round(optimal_rates, 0).astype(int),
        "% of Max": np.round(optimal_rates / max_oil * 100, 1),
        "Gas (Mscf/d)": np.round(optimal_rates * gor / 1000, 1),
        "Water (bbl/d)": np.round(optimal_rates * wc / (1 - wc), 0).astype(int),
        "Liquid (bbl/d)": np.round(optimal_rates / (1 - wc), 0).astype(int),
    })

    print("OPTIMAL PRODUCTION ALLOCATION")
    print("=" * 75)
    print(results.to_string(index=False))
    print("=" * 75)
    print(f"Total oil:    {total_oil:,.0f} STB/day")
    print(f"Total liquid: {sum(optimal_rates / (1-wc)):,.0f} bbl/day  (limit: {max_liquid_bpd:,})")
    print(f"Total gas:    {sum(optimal_rates * gor / 1e6):,.3f} MMscf/day (limit: {max_gas_mmscfd})")
    print(f"Total water:  {sum(optimal_rates * wc / (1-wc)):,.0f} bbl/day  (limit: {max_water_bpd:,})")

    # Compare to naive allocation (each well at max)
    naive_oil = sum(max_oil)
    print(f"\nNaive total (all at max): {naive_oil:,} STB/day")
    print(f"Optimized total:          {total_oil:,.0f} STB/day")
    if total_oil < naive_oil:
        print(f"Reduction due to constraints: {naive_oil - total_oil:,.0f} STB/day")
        print(f"But this allocation respects all facility limits.")
else:
    print(f"Optimization failed: {result.message}")

Gas Lift Optimization

Gas lift optimization across a field is a classic nonlinear allocation problem. Each well responds to gas injection with diminishing returns --- the first 500 Mscf/d of lift gas might add 300 barrels of oil, but the next 500 Mscf/d might only add 100. The total lift gas available is limited by compressor capacity. The question is: how to distribute the available gas across wells to maximize total oil production.

main.py

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

def gas_lift_response(q_base, gas_inj, a, b):
    """
    Gas lift performance curve.

    Returns the total oil rate as a function of gas injection rate.
    The response follows a diminishing-returns curve:
        q_oil = q_base + a * gas_inj / (b + gas_inj)

    Parameters
    ----------
    q_base : float
        Natural flow rate without gas lift (STB/day).
    gas_inj : float
        Gas injection rate (Mscf/day).
    a : float
        Maximum incremental oil gain (STB/day).
    b : float
        Gas rate at half-maximum response (Mscf/day).
    """
    return q_base + a * gas_inj / (b + gas_inj)


# Well gas lift response parameters
gl_wells = pd.DataFrame({
    "Well": ["GL-01", "GL-02", "GL-03", "GL-04"],
    "q_base": [400, 200, 600, 350],     # natural rate, STB/day
    "a": [800, 500, 400, 700],           # max incremental gain
    "b": [300, 200, 500, 250],           # half-response gas rate
})

# Total available lift gas
total_gas_mscfd = 2000  # Mscf/day from compressor

# Objective: maximize total oil (minimize negative)
def neg_total_oil(gas_alloc):
    total = 0
    for i, row in gl_wells.iterrows():
        total += gas_lift_response(row["q_base"], gas_alloc[i], row["a"], row["b"])
    return -total

# Constraints: total gas <= available
constraints = [{"type": "ineq", "fun": lambda x: total_gas_mscfd - sum(x)}]
bounds = [(0, 1500) for _ in range(len(gl_wells))]
x0 = [total_gas_mscfd / len(gl_wells)] * len(gl_wells)  # equal split initial guess

result = minimize(neg_total_oil, x0, method='SLSQP', bounds=bounds, constraints=constraints)

if result.success:
    opt_gas = result.x
    opt_oil = -neg_total_oil(opt_gas)

    # Compare to equal distribution
    equal_gas = [total_gas_mscfd / len(gl_wells)] * len(gl_wells)

    print("GAS LIFT OPTIMIZATION RESULTS")
    print("=" * 65)
    print(f"{'Well':<8} {'Base Rate':>10} {'Gas (Mscf/d)':>13} {'Oil Rate':>10} {'Gain':>8}")
    print("-" * 65)
    total_opt_oil = 0
    total_equal_oil = 0
    for i, row in gl_wells.iterrows():
        q_opt = gas_lift_response(row["q_base"], opt_gas[i], row["a"], row["b"])
        q_equal = gas_lift_response(row["q_base"], equal_gas[i], row["a"], row["b"])
        gain = q_opt - row["q_base"]
        total_opt_oil += q_opt
        total_equal_oil += q_equal
        print(f"{row['Well']:<8} {row['q_base']:>10,} {opt_gas[i]:>13,.0f} {q_opt:>10,.0f} {gain:>+8,.0f}")
    print("-" * 65)
    print(f"{'TOTAL':<8} {gl_wells['q_base'].sum():>10,} {sum(opt_gas):>13,.0f} {total_opt_oil:>10,.0f}")
    print(f"\nEqual distribution total: {total_equal_oil:,.0f} STB/day")
    print(f"Optimized total:         {total_opt_oil:,.0f} STB/day")
    print(f"Incremental gain:        {total_opt_oil - total_equal_oil:,.0f} STB/day")

    # Plot individual well responses
    fig, axes = plt.subplots(2, 2, figsize=(10, 8))
    gas_range = np.linspace(0, 1500, 200)

    for i, (ax, (_, row)) in enumerate(zip(axes.flat, gl_wells.iterrows())):
        oil_curve = [gas_lift_response(row["q_base"], g, row["a"], row["b"]) for g in gas_range]
        ax.plot(gas_range, oil_curve, 'k-', linewidth=2)
        ax.axvline(x=opt_gas[i], color='#2E86AB', linestyle='--', linewidth=1.5,
                   label=f'Optimal: {opt_gas[i]:.0f} Mscf/d')
        ax.axvline(x=equal_gas[i], color='#D7263D', linestyle=':', linewidth=1.5,
                   label=f'Equal: {equal_gas[i]:.0f} Mscf/d')
        ax.set_title(f'{row["Well"]}', fontsize=11)
        ax.set_xlabel('Gas Injection (Mscf/d)', fontsize=9)
        ax.set_ylabel('Oil Rate (STB/d)', fontsize=9)
        ax.legend(fontsize=8)
        ax.grid(True, alpha=0.3)

    plt.suptitle('Gas Lift Response Curves — Optimal vs. Equal Allocation', fontsize=13, y=1.01)
    plt.tight_layout()
    plt.show()

The optimizer allocates more gas to wells with steeper response curves (higher marginal oil gain per unit of gas) and less to wells that are already near their plateau. This is the principle of marginal analysis: distribute the resource to where it produces the greatest incremental return.

Production Surveillance

Optimization is not a one-time exercise. Reservoir conditions change, wells water out, equipment degrades. Production engineers need automated tools that continuously monitor well performance and flag anomalies.

main.py

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Generate synthetic production data for surveillance demo
np.random.seed(42)
days = 365
dates = pd.date_range("2025-01-01", periods=days, freq="D")

# Well with normal decline + sudden drop (possible problem)
base_rate = 1200 * np.exp(-0.0005 * np.arange(days))
noise = np.random.normal(0, 30, days)
normal_production = base_rate + noise

# Inject an anomaly: sudden drop at day 200, partial recovery
anomaly = normal_production.copy()
anomaly[200:230] = anomaly[200:230] * 0.4   # sudden 60% drop
anomaly[230:260] = anomaly[230:260] * 0.7   # partial recovery
anomaly[260:] = anomaly[260:] * 0.85         # new lower baseline

prod_data = pd.DataFrame({
    "Date": dates,
    "Oil Rate (STB/d)": np.round(anomaly, 0),
})

# Rolling statistics for anomaly detection
window = 14  # 14-day rolling window
prod_data["Rolling Mean"] = prod_data["Oil Rate (STB/d)"].rolling(window).mean()
prod_data["Rolling Std"] = prod_data["Oil Rate (STB/d)"].rolling(window).std()
prod_data["Upper Bound"] = prod_data["Rolling Mean"] + 2 * prod_data["Rolling Std"]
prod_data["Lower Bound"] = prod_data["Rolling Mean"] - 2 * prod_data["Rolling Std"]

# Flag anomalies (rate outside 2-sigma bounds)
prod_data["Anomaly"] = (
    (prod_data["Oil Rate (STB/d)"] < prod_data["Lower Bound"]) |
    (prod_data["Oil Rate (STB/d)"] > prod_data["Upper Bound"])
)

anomaly_count = prod_data["Anomaly"].sum()
anomaly_dates = prod_data[prod_data["Anomaly"]]["Date"]

print("PRODUCTION SURVEILLANCE REPORT")
print("=" * 50)
print(f"  Monitoring period: {dates[0].date()} to {dates[-1].date()}")
print(f"  Total days:        {days}")
print(f"  Anomalies flagged: {anomaly_count}")
print(f"  First anomaly:     {anomaly_dates.iloc[0].date() if len(anomaly_dates) > 0 else 'None'}")
print()

# Plot
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(prod_data["Date"], prod_data["Oil Rate (STB/d)"], 'k-', linewidth=0.8, alpha=0.7, label='Daily Rate')
ax.plot(prod_data["Date"], prod_data["Rolling Mean"], 'b-', linewidth=1.5, label=f'{window}-day Mean')
ax.fill_between(prod_data["Date"], prod_data["Lower Bound"], prod_data["Upper Bound"],
                alpha=0.15, color='blue', label='2σ Bounds')
anomaly_pts = prod_data[prod_data["Anomaly"]]
ax.scatter(anomaly_pts["Date"], anomaly_pts["Oil Rate (STB/d)"],
           color='red', s=15, zorder=5, label='Anomaly')
ax.set_xlabel('Date', fontsize=11)
ax.set_ylabel('Oil Rate (STB/day)', fontsize=11)
ax.set_title('Well A-01 — Production Surveillance', fontsize=12)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Exercises

fitness_center

Exercise 13.1Practice

-- Basic Rate Allocation

Three wells with maximum rates of 1,000, 800, and 1,500 STB/day share a pipeline with 2,500 STB/day capacity. Each well has operating costs of 10,10, ...

arrow_forward