How to Predict Oilfield Production Using Grey Models and Python

Oilfield Prediction

Problem

The S reservoir of a Chinese oilfield has monthly production data from February 1994 to February 1995. The task is to build a predictive model for monthly oil production, water production, water injection, and bottom‑hole pressure for March 1995.

Modeling

Let x1, x2, x3, x4 denote the monthly oil production, water production, water injection, and bottom‑hole pressure, respectively. A GM(1,1) model is built for oil production, a GM(1,2) model linking water production to oil production, a GM(1,1) model for water injection, and a GM(1,3) model relating pressure to oil and water injection. The equations are rewritten in matrix form and the coefficients are identified using least‑squares.

Solution

The system of differential equations is solved numerically with Python (odeint). Using the February 1994 values as initial conditions, the model generates predictions for March 1995.

Verification

Simulation results match the observed values reasonably well, except for a relatively large error in one water‑injection prediction. The predicted March 1995 values are presented in the program output.

Code

<code>import numpy as np
from scipy.integrate import odeint
a=np.loadtxt("Pdata15_3.txt")  # load the last four columns
n=a.shape[0]  # number of observations
x10=a[:,0]; x20=a[:,1]; x30=a[:,2]; x40=a[:,3]
x11=np.cumsum(x10); x21=np.cumsum(x20)
x31=np.cumsum(x30); x41=np.cumsum(x40)

z1=(x11[:-1]+x11[1:])/2.; z2=(x21[:-1]+x21[1:])/2.
z3=(x31[:-1]+x31[1:])/2.; z4=(x41[:-1]+x41[1:])/2.
B1=np.c_[z1,np.ones((n-1,1))]
u1=np.linalg.pinv(B1).dot(x10[1:]); print(u1)
B2=np.c_[z1,z2]
u2=np.linalg.pinv(B2).dot(x20[1:]); print(u2)
B3=np.c_[z3,np.ones((n-1,1))]
u3=np.linalg.pinv(B3).dot(x30[1:]); print(u3)
B4=np.c_[z1,z3,z4]
u4=np.linalg.pinv(B4).dot(x40[1:]); print(u4)

def Pfun(x,t):
    x1, x2, x3, x4 = x
    return np.array([u1[0]*x1+u1[1],
                     u2[0]*x1+u2[1]*x2,
                     u3[0]*x3+u3[1],
                     u4[0]*x1+u4[1]*x3+u4[2]*x4])

t=np.arange(0,14)
X0=np.array([7.1230,0.7960,13.1080,27.475])
s1=odeint(Pfun, X0, t); s2=np.diff(s1,axis=0)
xh=np.vstack([X0,s2])
cha=a-xh[:-1,:]  # residuals
delta=np.abs(cha/a)  # relative errors
maxd=delta.max(0)  # maximum relative error for each indicator
pre=xh[-1,:]; print("Maximum relative error:",maxd,"\nPrediction:",pre)
</code>

Reference

Si Shoukuai, Sun Xijing, Python Mathematics Experiments and Modeling