数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

油田预测

问题

我国某油田 S 油藏自 1994 年 2 月至 1995 年 2 月的开发动态 数据见下表, 试建立预测模型, 并预测月产油量、月产水量、月注水量和 底层压力 1995 年 3 月的取值。

序号 时间 月产油量/万吨 月产水量/万吨 月注水量/万吨 底层压力/MPa
1 94.02 7.123 0.796 13.108 27.475
2 94.03 7.994 0.832 12.334 27.473
3 94.04 8.272 0.917 12.216 27.490
4 94.05 7.960 0.976 12.201 27.500
5 94.06 7.147 1.075 12.132 27.510
6 94.07 7.092 1.121 11.990 27.542
7 94.08 6.858 1.281 11.926 27.536
8 94.09 5.804 1.350 10.478 27.550
9 94.10 6.433 1.410 9.176 27.567
10 94.11 6.354 1.432 11.368 27.584
11 94.12 6.254 1.507 12.764 27.600
12 95.01 5.197 1.559 11.143 27.602
13 95.02 5.654 1.611 10.737 27.630

建模

用 $x_{1}, x_{2}, x_{3}, x_{4}$ 分别表示月产油量、月产水量、月注水量和底层压力, 它们的观测值分别记作 $x_{i}^{(0)}=\left(x_{i}^{(0)}(1), x_{i}^{(0)}(2), \cdots, x_{i}^{(0)}(13)\right)(i=1,2,3,4)$ 。

对月产油量 $x_{1}$, 建立 GM(1,1)模型; 月产水量 $x_{2}$ 受月产油量及油藏水驱 规律的控制, 因此对 $x_{2}$ 建立与 $x_{1}$ 相关的 GM(1,2)模型; 对月注水量 $x_{3}$, 建立 $\mathrm{GM}(1,1)$ 模型; 地层压力 $x_{4}$ 与采注关系密切, 是地层能量的综合反映, 因此对 $x_{4}$ 建立与 $x_{1} 、 x_{3}$ 相关的 GM(1,3)模型。综上所述, 系统的 4 个开发指标所满足的状态方程为

$$ \left\{\begin{array}{l} \frac{d x_{1}^{(1)}(t)}{d t}=a_{11} x_{1}^{(1)}(t)+b_{1} \\ \frac{d x_{2}^{(1)}(t)}{d t}=a_{21} x_{1}^{(1)}(t)+a_{22} x_{2}^{(1)}(t) \\ \frac{d x_{3}^{(1)}(t)}{d t}=a_{33} x_{3}^{(1)}(t)+b_{3} \\ \frac{d x_{4}^{(1)}(t)}{d t}=a_{41} x_{1}^{(1)}(t)+a_{43} x_{3}^{(1)}(t)+a_{44} x_{4}^{(1)}(t) \end{array}\right. $$

将式子改写为矩阵方程 $$ \frac{d X^{(1)}(t)}{d t}=A X^{(1)}(t)+\tilde{B} $$

式中, $$ \begin{aligned} \frac{d X^{(1)}(t)}{d t} &=\left[\frac{d x_{1}^{(1)}(t)}{d t}, \frac{d x_{2}^{(1)}(t)}{d t}, \frac{d x_{3}^{(1)}(t)}{d t}, \frac{d x_{4}^{(1)}(t)}{d t}\right]^{T}, \\ X^{(1)}(t)=& {\left[\begin{array}{l} x_{1}^{(1)}(t) \\ x_{2}^{(1)}(t) \\ x_{3}^{(1)}(t) \\ x_{4}^{(1)}(t) \end{array}\right], \quad A=\left[\begin{array}{cccc} a_{11} & 0 & 0 & 0 \\ a_{21} & a_{22} & 0 & 0 \\ 0 & 0 & a_{33} & 0 \\ a_{41} & 0 & a_{43} & a_{44} \end{array}\right], \quad \tilde{B}=\left[\begin{array}{c} b_{1} \\ 0 \\ b_{3} \\ 0 \end{array}\right] } \end{aligned} $$ 利用上面介绍的最小二乘法, 依次辨识每个方程中的系数 $$ a_{11}, b_{1} ; a_{21}, a_{22} ; a_{33}, b_{3} ; a_{41}, a_{43}, a_{44} \circ $$

计算 4 茨得到 $$ A=\left[\begin{array}{cccc} -0.0378 & 0 & 0 & 0 \\ 0.0568 & -0.2232 & 0 & 0 \\ 0 & 0 & -0.0114 & 0 \\ -0.6179 & 0 & 5.6950 & -2.1926 \end{array}\right], \quad \tilde{B}=\left[\begin{array}{c} 8.6667 \\ 0 \\ 12.4938 \\ 0 \end{array}\right] $$

以表中 1994 年 2 月的月产油量、月产水量、月注水量、地层压力的数据为初值, 得如下初值问题: $$ \left\{\begin{array}{l} \frac{d X^{(1)}(t)}{d t}=A X^{(1)}(t)+\tilde{B}, \\ X^{(1)}(1)=[7.123,0.796,13.108,27.475] . \end{array}\right. $$

求解

利用 Python 软件求微分方程组 的数值解, 可以得到 $X^{(1)}(k)$ $(k=2,3, \cdots, 14)$ 的数值解。

检验

从程序运行结果可以看出, 模拟值与实际值符合较好, 但有一个注水值 的预测值的相对误差较大。1995 年 3 月 $x_{1}, x_{2}, x_{3}, x_{4}$ 的预测值分别为 $5.2347$, $$ 1.4599,10.7070,26.5756 \text { 。 } $$

代码

import numpy as np
from scipy.integrate import odeint
a=np.loadtxt("Pdata15_3.txt")  #加载表中的后4列数据
n=a.shape[0]  #观测数据的个数
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,:]  #计算残差

delta=np.abs(cha/a)  #计算相对误差
maxd=delta.max(0)  #计算每个指标的最大相对误差
pre=xh[-1,:]; print("最大相对误差：",maxd,"\n预测值为：",pre)