Well control optimization using a two-step surrogate treatment

TL;DR

This paper introduces a two-step surrogate treatment method for well control optimization in waterflood oil production, significantly reducing computational costs while maintaining near-optimal results, suitable for iterative field development planning.

Contribution

The paper presents a novel two-step surrogate treatment that efficiently approximates well control optimization, enabling faster computations with minimal accuracy loss compared to traditional methods.

Findings

Achieves 5x or more speedup over formal optimization methods.

Provides results in reasonable agreement with MADS and adjoint-gradient.

Applicable for 2D and 3D models with different fluid mobilities.

Abstract

Large numbers of flow simulations are typically required for determining optimal well settings. These simulations are often computationally demanding, which poses challenges for the optimizations. In this paper we present a new two-step surrogate treatment (ST) that reduces the computational expense associated with well control optimization. The method is applicable for oil production via waterflood, with well rates optimized at a single control period. The two-step ST entails two separate optimizations, which can both be performed very efficiently. In the first optimization, optimal well-rate ratios (i.e., the fraction of total injection or production associated with each well) are determined such that a measure of velocity variability over the field is minimized, leading to more uniform sweep. In the second step, overall injection and production rates are determined. The flow physics in the first step is highly simplified, while the actual physical system is simulated in the second step. Near-globally-optimal results can be determined in both cases, as the first optimization is posed as a QP problem, and the second step entails just a single optimization variable. Under full parallelization, the overall elapsed time for the ST corresponds to the runtime for 1-2 full-order simulations. Results are presented for multiple well configurations, for 2D and 3D channelized models, and comparisons with formal optimization procedures (MADS and adjoint-gradient) are conducted. Three different fluid mobility ratios (M = 1, 3 and 5) are considered. Optimization results demonstrate that the two-step ST provides results in reasonable agreement with those from MADS and adjoint-gradient, with speedups of 5x or more. We also show that the ST is applicable in the inner-loop in field development optimization, where it will be especially useful since many different well configurations must be evaluated.

Tables

Algorithm 1. Well-rate ratios and field rate optimization
1	Solve first-step optimization problem to provide ${\text{f}}^{\ast }$ (Matlab/CVX) (Eq. 16)
2	Solve second-step optimization problem to provide $Q^{\ast }$ (Eq. 17):
3	for number of computing nodes (or licences) available do
4	Run full-physics simulation with different field rates (AD-GPRS)
5	Treat BHP constraints as described in Sect. 3.2
6	end for
7	Choose field rate associated with largest NPV

Table 2. Table 1 : Simulation and optimization parameters for 100 × 100 100 100 100\times 100 channelized model

Grid dimensions	100 $\times$ 100 $\times$ 1
Grid cell dimensions	50 ft $\times$ 50 ft $\times$ 100 ft
Initial pressure $p_i$, at datum	6000 psi at 6000 ft
${\mu }_o$	1.0, 3.0 or 5.0 cp
${\mu }_w$	1.0 cp
$B_o$ and $B_w$ at $p_i$	1.075 RB/STB
Oil compressibility at $p_i$	10-5 psi-1
Water compressibility at $p_i$	10-5 psi-1
Rock compressibility at $p_i$	10-9 psi-1
Porosity	0.25
Simulation period	7300 days
$p_o$, $c_{pw}$ and $c_{iw}$	$60, $5 and $5/STB
$c_w$	$20 MM
Discount rate	0$\%$

Table 3. Table 2 : NPVs for base cases ( M = 3 𝑀 3 M=3 ) for 2D reservoir model

Well	NPV (BN $)
scenario
	PVI = 0.25	0.75	1.25
$1^1$	0.54	0.48	0.15
$2^2$	0.81	0.89	0.62
$3^3$	0.49	1.22	1.41
${}^1{Q}_{1,M=3}^{\ast }=0.49\text{PVI}$
${}^2{Q}_{2,M=3}^{\ast }=0.41\text{PVI}$
${}^3{Q}_{3,M=3}^{\ast }=0.65\text{PVI}$

Table 4. Table 3 : Optimal NPVs using different methods (2D reservoir model)

Well	Mobility	NPV (BN $)
scenario	ratio			Adjoint-
		ST	${\text{MADS}}^1$	${\text{gradient}}^1$
	1	1.23	1.24	1.24
1	3	0.62	0.62	0.62
	5	0.47	0.47	0.47
	1	2.03	2.02	2.07
2	3	1.03	1.09	1.02
	5	0.84	0.89	0.89
	1	2.55	2.64	2.55
3	3	1.47	1.47	1.47
	5	1.28	1.29	1.29
${}^1$Median values

Table 5. Table 4 : Minimum and maximum NPVs for MADS and adjoint-gradients (2D reservoir model)

Well	Mobility	NPV (BN $)
scenario	ratio				Adjoint-
		ST	MADS		gradient
			min	max	min	max
	1	1.23	1.23	1.24	1.24	1.24
1	3	0.62	0.61	0.62	0.62	0.62
	5	0.47	0.46	0.47	0.46	0.47
	1	2.03	2.02	2.02	2.07	2.07
2	3	1.03	1.09	1.09	1.02	1.03
	5	0.84	0.89	0.89	0.87	0.89
	1	2.55	2.64	2.64	2.55	2.55
3	3	1.47	1.47	1.47	1.43	1.47
	5	1.28	1.28	1.29	1.11	1.29

Table 6. Table 5 : Elapsed time (fully parallelized for ST and MADS) for different methods for 2D reservoir model

Well	Mobility	Elapsed time (s)
scenario	ratio			Adjoint-
		ST	${\text{MADS}}^1$	${\text{gradient}}^1$
	1	52	389	1652
1	3	52	385	1952
	5	51	389	3665
	1	60	361	1138
2	3	63	458	1165
	5	62	456	2538
	1	67	503	876
3	3	68	503	2750
	5	66	502	3690
${}^1$Median values

Table 7. Table 6 : Simulation and optimization parameters for 40 × 30 × 9 40 30 9 40\times 30\times 9 channelized model

Grid dimensions	40 $\times$ 30 $\times$ 9
Grid cell dimensions	100 ft $\times$ 100 ft $\times$ 15 ft
Initial pressure $p_i$, at datum	6000 psi at 6000 ft
${\mu }_o$	1.0, 3.0 or 5.0 cp
${\mu }_w$ (constant)	1.0 cp
${\rho }_o$ and ${\rho }_w$	53.10 and 64.79 lbm/ft3
$B_o$ and $B_w$ at $p_i$	1.075 RB/STB
Oil compressibility at $p_i$	10-5 psi-1
Water compressibility at $p_i$	10-5 psi-1
Rock compressibility at $p_i$	10-9 psi-1
Field porosity	0.25
Simulation period	14,600 days
$p_o$, $c_{pw}$ and $c_{iw}$	$60, $5 and $5/STB
$c_w$	$20 MM
Discount rate	0% or 5%

Table 8. Table 7 : NPVs for base cases ( M = 3 𝑀 3 M=3 ) for 3D reservoir model with a single control period

Well	NPV (BN $)
scenario
	PVI = 1.0	1.5	2.0
$1^1$	3.71	4.72	5.22
$2^2$	3.69	5.14	6.61
${}^1{Q}_{1,M=3}^{\ast }=1.44\text{PVI}$
${}^2{Q}_{2,M=3}^{\ast }=1.84\text{PVI}$

Table 9. Table 8 : Optimal NPVs for ST and MADS for 3D reservoir model with a single control period

Well	Mobility	NPV (BN $)
scenario	ratio		MADS
		ST	med	min	max
	1	17.64	17.44	17.35	17.53
1	3	5.45	5.52	5.48	5.60
	5	3.85	4.04	4.04	4.04
	1	22.70	22.88	22.88	22.88
2	3	7.32	7.56	7.56	7.57
	5	4.65	4.80	4.80	4.80

Table 10. Table 9 : Elapsed time (fully parallelized) for ST and MADS for a single control period for 3D reservoir model

Well	Mobility	Elapsed time (s)
scenario	ratio
		ST	${\text{MADS}}^1$
	1	151	741
1	3	146	736
	5	147	738
	1	132	741
2	3	136	735
	5	137	734
${}^1$Median values

Table 11. Table 10 : Optimal NPVs for 3D reservoir model with multiple control periods (well scenario 1).

		NPV (BN $)
Mobility	d	ST	ST &	MADS
ratio	(%)		MADS	med	min	max
3	0	5.45	5.63	5.46	5.36	5.55
5	5	1.71	1.76	1.76	1.46	1.76

Table 12. Table 11 : Elapsed time (fully parallelized) for 3D reservoir model with multiple control periods (well scenario 1).

	Elapsed time (s)
Mobility	ST	ST &	MADS
ratio		MADS	med	min	max
3	136	971	1355	1228	1377
5	137	1211	1515	1493	1655

Equations

$$\nabla \cdot \left({\lambda }_t\left(S_w\right)\text{k}\nabla p\right)=-{\tilde{q}}_t^w+\varphi c_t\frac{\partial p}{\partial t},$$

$$\varphi \frac{\partial S_w}{\partial t}+\nabla \cdot \left(f\left(S_w\right){\text{u}}_t\right)={\tilde{q}}_w^w.$$

$$\displaystyle\begin{array}[]{rrclcl}\displaystyle\min_{\textbf{x}\in\mathbb{X},\textbf{u}\in\mathbb{U},\textbf{z}\in\mathbb{Z}}&{J(\textbf{x},\textbf{u},\textbf{z})},\ \ \textrm{subject to}\begin{cases}\textbf{g}(\textbf{p},\textbf{x},\textbf{u},\textbf{z})=\textbf{0},\\ \textbf{c}(\textbf{p},\textbf{x},\textbf{u},\textbf{z})\leq\textbf{0}\end{cases}\end{array}$$

$$\text{NPV(x,u,z)}=\sum _{k=1}^{n_p}\sum _{s=1}^{n_s}\frac{\Delta t_s\left(p_o{q}_{k,s}^o(\text{x},\text{u})-c_{pw}{q}_{k,s}^{pw}(\text{x},\text{u})\right)}{(1+d{)}^{\frac{t_s}{365}}}-$$

$$\sum _{k=1}^{n_i}\sum _{s=1}^{n_s}\frac{\Delta t_s{c}_{iw}{q}_{k,s}^{iw}(\text{x},\text{u})}{(1+d{)}^{\frac{t_s}{365}}}-\sum _{k=1}^{n_w}\frac{|z_k|c_w}{(1+d{)}^{\frac{t_k}{365}}}.$$

$$\begin{array}[]{rrclcl}\displaystyle\max_{\textbf{u}\in\mathbb{U}}&{\textnormal{NPV}(\textbf{u})},\ \ \textrm{subject to}\begin{cases}\textbf{g}(\textbf{p},\textbf{u})=\textbf{0},\\ \textbf{c}(\textbf{p},\textbf{u})\leq\textbf{0}\end{cases}\end{array}$$

$$\text{qˆ}=\text{f}Q,$$

$${\left[{{\text{f}}^{\ast }}^{\text{T}},Q^{\ast }\right]}^{\text{T}}=\underset{\text{f}\in F,Q\in \mathbb{R}}{\mathrm{argmin}}J(\text{f},Q),$$

$${\text{f}}^{\ast }=\underset{\text{f}\in F}{\mathrm{argmin}}{J}_a(\text{f}),$$

$$Q^{\ast }=\underset{Q\in \mathbb{R}}{\mathrm{argmin}}{J}_b({\text{f}}^{\ast },Q),$$

$$\nabla \cdot \left(\frac{\text{k}}{\mu }\nabla p\right)=-{\tilde{q}}_t^w+\varphi c\frac{\partial p}{\partial t},$$

$$\text{u}=-\frac{\text{k}}{\mu }\nabla p.$$

$$(u_x{)}_{i\pm \frac{1}{2},j,k}=\frac{(\tilde{q}{)}_{i\pm \frac{1}{2},j,k}}{A}.$$

$$(u_x{)}_{i,j,k}=\frac{(u_x{)}_{i-\frac{1}{2},j,k}+(u_x{)}_{i+\frac{1}{2},j,k}}{2}.$$

$${\text{uˆ}}_d={\text{V}}_d\text{qˆ},$$

$$\underset{\text{f}\in F}{\min }{\text{uˆ}}_x^{\text{T}}{\text{W}}_x{\text{uˆ}}_x+{\text{uˆ}}_y^{\text{T}}{\text{W}}_y{\text{uˆ}}_y+{\text{uˆ}}_z^{\text{T}}{\text{W}}_z{\text{uˆ}}_z,$$

$$\text{subject to}\sum _{k=1}^{n_p}{f}_k^p=\sum _{k=1}^{n_i}{f}_k^i=1,f_k^p,f_k^i\ge 0,f\in \mathbb{R},$$

$$\underset{Q\in \mathbb{R}}{\max }\text{NPV}({\text{f}}^{\ast },Q),\text{subject to}Q>0.$$

$$(f^{\prime }{)}_k^p=\frac{f_k^p}{\sum _{k=1}^{n_{pr}}{f}_k^p},$$

Full text

††journal: my journal

1 Introduction

Well control optimization entails defining the set of controls (well rates or bottom-hole pressures) that maximize a predefined cost function, such as net present value or total oil production. Optimization computations generally rely on a reservoir simulator for the evaluation of proposed well control strategies. Since the flow models are typically time dependent and nonlinear, the solution of the resulting nonconvex optimization problem may be computationally demanding. This is even more of an issue when well locations and controls are optimized simultaneously. In such cases, the objective function evaluation for each set of well locations requires the determination of optimal well controls.

In this study, we introduce a new two-step surrogate treatment (ST) for well control optimization for oil fields producing under waterflood. In the first step of the ST, we determine the well-rate ratios (i.e., the fraction of total injection or production associated with each well), for an appropriately defined unit-mobility-ratio problem. In the second step of the ST, we fix the well-rate ratios found in the first step, and then solve a full-physics optimization problem in one variable to determine the overall injection/production rate that maximizes net present value (NPV). Although the problems in the ST are idealized, they are posed such that we can efficiently find near-global optima in both steps (though this does not mean we determine the global optimum for the actual problem).

Well control optimization problems have been addressed using both adjoint-gradient and derivative-free methods. Adjoint-gradient-based methods, which include those described in [1, 2, 3], are very efficient though they require access to reservoir simulation source code. Derivative-free optimization (DFO) methods represent a useful alternative when gradients are not available. Both stochastic and pattern search DFO methods have been applied in well control optimization problems. Stochastic methods used in this setting include genetic algorithms (GAs) and particle swarm optimization (PSO), which are both population-based evolutionary methods. These procedures were shown to be applicable for both well control [4, 5, 6, 7] and well location optimization [8, 9]. Pattern-search or stencil-based optimization techniques, such as mesh adaptive direct search (MADS) [10], have also been used for well control optimization [11, 12, 13, 14]. In this work both adjoint-gradient-based and MADS procedures will be compared to our new two-stage surrogate treatment.

The use of proxy or surrogate models is an effective way to reduce the cost of optimization in cases when the underlying full-physics model is expensive to simulate. Surrogate treatments have been applied for both well control and well placement optimization problems. For example, in an early study involving field development optimization, Rosenwald and Green [15] used single-phase flow models, for which linear programming techniques are applicable. Subsequent studies that used surrogate models in the form of statistical proxies, for well placement optimization, include those of Guyaguler et al. [16], Yeten et al. [17] and Farmer et al. [18]. Reduced-physics surrogate models were used by Wilson and Durlofsky [19] for optimizing the development of shale gas reservoirs [19], while Aliyev and Durlofsky [20] applied a sequence of upscaled models as surrogates for joint (well location and control) optimization in waterflooding problems. Coarse models were also used as surrogates in well placement optimization by Bukhamsin [21].

Surrogate treatments for well control optimization, based on reduced-order modeling with proper orthogonal decomposition, have been developed by a number of researchers; see, e.g., [22, 23, 24], along with the review by Jansen and Durlofsky [25]. Reduced-physics models based on streamline methods were shown to provide useful surrogates in water injection optimization [26]. Well control optimization applications using streamlines have also been described by Park and Datta-Gupta [27], who used 3D streamlines to generate flood efficiency maps, and by Wen et al. [28], who applied both streamline and time-of-flight (TOF) computations for the optimization of mature fields under waterflood. Flow diagnostics tools [29, 30], typically constructed based on reduced-physics simulations, represent another effective set of surrogate treatments that can be used for optimization. For example, M$\o$yner et al. [31] performed well control optimization for systems involving two- and three-phase flow. Lie et al. [30] used quantities derived from TOF and multiscale methods as fast proxies for well control optimization under waterflood. We note finally that He et al. [32] and Chen et al. [33] have demonstrated that surrogate models can also be applied to accelerate history matching and uncertainty quantification computations.

Our two-step ST somewhat resembles the procedure developed by Rodríguez Torrado et al. [34] and Embid Droz et al. [35] for setting well controls. Their methodology also relies on the optimization of a single overall field rate variable (e.g., field oil production rate), and was used in combination with well scheduling and location optimization procedures. Well controls with this approach are determined through estimation of the field recovery factor, while the field rate variable is obtained using derivative-free methods. A key difference between our two-step ST and the procedure for optimizing the well controls developed by Rodríguez Torrado et al. [34] is in the (heuristic) scheme used to determine the well-rate ratios. In our approach, a simplified flow problem is considered for the determination of the well-rate ratios, while their procedure uses static reservoir properties. Detailed comparisons would, however, need to be performed to clearly establish the relative advantages of the two methods.

In this paper, we develop and apply a new two-step surrogate treatment (ST) to enable efficient well control optimization for subsurface flow problems. Although field development optimization is not the focus of this paper, a key applications for our ST is the joint optimization of well locations and controls, where the ST can be used as a fast optimization method for the (inner-loop) well control optimization. The computational savings achieved by ST derive from solving two approximate subproblems rather than the actual full-physics problem. In the first subproblem we simplify the physics, and in the second subproblem we consider the actual physics but greatly simplify the optimization problem. Optimum solutions can be obtained for both subproblems at a cost of 1–2 full-order runs under full parallelization. The method provides well controls for only a single control period, but this solution can then be used as an initial guess for cases involving multiple control periods. Geological uncertainty is not considered, though this could be addressed by optimizing over multiple realizations.

This paper is organized as follows. In Sect. 2 we present the basic oil-water flow equations and the general well control optimization problem. In Sect. 3, the two-step surrogate treatment is described in detail. Heuristic procedures for handling BHP constraints are discussed, and the general workflow is explained. In Sect. 4, we present ST well control optimization results for 40 different (waterflood) well scenarios in a 2D model. Well control optimization for additional mobility ratios, and comparisons with optimization results using adjoint-gradient and MADS procedures, are also provided. Results for a 3D case are then presented. Finally, in Sect. 5, we summarize this work and provide suggestions for future research in this area.

2 Problem statement

In this section, we present the equations governing reservoir flow, and we provide a formal statement of the general optimization problem.

2.1 Oil-water flow equations

For two-phase oil-water flow, the governing equations, which derive from Darcy’s law for each phase and mass conservation statements, can be expressed in terms of the so-called pressure and saturation equations. In the absence of capillary pressure and gravity effects, these equations are as follows:

[TABLE]

[TABLE]

Here subscripts $o$ and $w$ refer to oil and water phases, $\lambda_{t}$ is the total mobility, with $\lambda_{t}=\lambda_{o}+\lambda_{w}$, where $\lambda_{o}=k_{ro}/\mu_{o}$ and $\lambda_{w}=k_{rw}/\mu_{w}$, $k_{ro}$ and $k_{rw}$ are relative permeability to the oil and water phases respectively, and $\mu_{o}$ and $\mu_{w}$ are phase viscosities. Additional variables are water saturation $S_{w}$, absolute permeability tensor k, pressure $p$, total source/sink $\tilde{q}^{w}_{t}$ (the superscript $w$ here refers to wells), porosity $\phi$, total compressibility $c_{t}$ (which is a weighted combination of the fluid and rock compressibilities), and time $t$. The fractional flow of water $f(S_{w})$ is defined as $f(S_{w})=\lambda_{w}/\lambda_{t}$, the total velocity as $\textbf{u}_{t}=\textbf{u}_{o}+\textbf{u}_{w}$ ($\textbf{u}_{o}$ and $\textbf{u}_{w}$ are the Darcy velocities for oil and water), and $\tilde{q}^{w}_{w}$ is the water source/sink term. We additionally have the saturation constraint $S_{o}+S_{w}=1$, which completes the problem statement.

In the discretized system, the source terms ($\tilde{q}^{w}_{w}$ and $\tilde{q}^{w}_{t}$) are represented using the Peaceman well model. These terms are nonzero only in grid blocks containing wells that are open to flow. The well model relates flow rate to bottom-hole pressure (BHP), and either well rates or BHPs can be specified as the control variables. In this study we consider vertical wells, completed over the entire reservoir thickness, with flow rates as the control parameter. Stanford’s Automatic Differentiation-based General Purpose Research Simulator (AD-GPRS) [36] is used for all flow simulations. Although gravitational terms are not shown in Eqs. 1 and 2, in order to simplify the presentation, these effects are included in the 3D simulations performed in Sect. 4.

2.2 General optimization problem

The general field development optimization problem involves determination of the well types, locations and controls, with the goal of minimizing a cost function $J$. Following [37], the optimization problem can be stated as follows:

[TABLE]

The vectors x and u indicate integer (grid-block based) well location variables and continuous well control variables, respectively, while z are categorical variables, which indicate whether the well is an injector ($z_{k}=-1$), a producer ($z_{k}=1$), or not drilled at all ($z_{k}=0$). The well location variables can also be treated as real-valued, and this may be preferable in cases where wells are not centered in grid blocks (e.g., with deviated wells). Here $\textbf{g}=\textbf{0}$ denotes the flow simulation equations, p represents the solution unknowns, which in our system are the pressure and saturation in every grid block, and c defines any nonlinear constraints. The spaces $\mathbb{X}$ and $\mathbb{U}$ are defined to include bound constraints, which can be expressed as $\textbf{x}_{l}\leq\textbf{x}\leq\textbf{x}_{u}$ and $\textbf{u}_{l}\leq\textbf{u}\leq\textbf{u}_{u}$, where subscripts $l$ and $u$ denote lower and upper bounds.

In this work the objective is to maximize net present value (NPV); i.e., we set $J=-{\textrm{NPV}}$, with NPV given by:

[TABLE]

Here $n_{i}$ is the number of injection wells, $n_{p}$ is the number of production wells, $n_{w}=n_{i}+n_{p}$ is the total number of wells, $n_{s}$ is the number of simulation time steps, $t_{s}$ and $\Delta t_{s}$ are the time and time step size at time step $s$, and $d$ is the annual discount rate. The price of oil, the cost for handling produced water, and the cost of injected water are, respectively, $p_{o}$, $c_{pw}$ and $c_{iw}$. The oil and water production rates and the water injection rate, for well $k$ at time step $s$, are denoted ${q}^{o}_{k,s}$, ${q}^{pw}_{k,s}$ and ${q}^{iw}_{k,s}$. The variable $t_{k}$ represents the time at which well $k$ is drilled, and the per-well drilling cost is denoted by $c_{w}$. For a particular well control optimization problem, the last term in Eq. 5 does not vary since the well configuration, and thus the values of $|z_{k}|$, are fixed. We nonetheless include the well costs in the computation in order to render the resulting NPVs more indicative of those for an actual project.

In the optimizations performed in this study, the well types and locations are fixed (though in one set of examples we extract the well configurations from a field development optimization run). Then, with $J=-$NPV(u), Eq. 4 can be expressed as:

[TABLE]

The space $\mathbb{U}$ again includes the bound constraints for the continuous well control variables, and c specifies any nonlinear constraints.

We note that, given the nature of the optimization problem, fixing the well locations and types represents a considerable simplification. This specification eliminates the integer well placement variables (these variables are integers because well locations are represented in terms of discrete grid blocks), as well as the categorical well type variables. If these variables are included in the optimization, the overall field development optimization problem (Eq. 4) is a mixed integer nonlinear programming (MINLP) problem, which is much more difficult to solve than the real-variable problem defined in Eq. 6.

The MINLP field development problem has been addressed using both sequential and joint optimization procedures. In the sequential approach, the well types and locations are optimized first, with well controls, or a well control strategy, specified. Once the optimized well placement configuration is determined, the well controls are then optimized for this configuration. However, since the optimal well locations depend on how the wells are operated, the sequential approach will in general provide a suboptimal solution. In the joint optimization approach, well types, locations, and controls are optimized together. This can be accomplished either in a single loop [38], in which all variables are considered simultaneously, or in a nested fashion [39]. In the nested approach, well types and locations are defined in the outer loop, and well controls are optimized in the inner-loop. The ST developed in this study can be used directly in the inner-loop in nested (joint) field development optimization procedures. In this more challenging application, our two-step ST may provide significant computational speedup.

It is very useful for the surrogate model to be constructed such that the global optimum of the resulting optimization problem can be readily found. This will be the case, for example, if the surrogate treatment leads to a linear programming (LP) or a quadratic programming (QP) problem. For general LP or QP problems, a variety of global optimization algorithms are available. These include the Simplex algorithm and its extensions, and interior-point and conjugate-gradient methods [40]. As we will see, the two-step ST developed in this work can indeed be formulated as two supbroblems that are amenable to global optimization. We reiterate that this does not mean we find the global optimum for the original full-physics problem, but rather that we can find the global optima for the subproblems considered in the ST.

3 Optimization methodology

In this section we present our two-step surrogate treatment for well control optimization of oil fields under waterflood. Both ST steps are described in detail, and the handling of BHP constraints is also discussed. This treatment considers only a single control period. Although our development here is for production via waterflood, ST is also expected to be applicable for other displacement processes, such as gas or water-alternating-gas injection, though it may be more approximate for these operations.

3.1 Two-step ST

We denote the optimization variables used in this problem as $\textbf{f}\in\mathbb{R}^{n_{w}\times 1}$, where $n_{w}$ is the number of wells, and $Q\in\mathbb{R}$. Here f represents the vector of well-rate ratios, or the fraction of the field rate that is allocated to each well, and Q specifies the total field injection and production rates. At this point these total rates are the same since we have unit voidage replacement ratio, VRR = 1. With these definitions, the actual well rates can be obtained from the following relationship between f and Q:

[TABLE]

where $\textbf{\^{q}}=\left[(\textbf{\^{q}\,\textsuperscript{\it{p}}})^{\textnormal{T}},(\textbf{\^{q}\,\textsuperscript{\it{i}}})^{\textnormal{T}}\right]^{\textnormal{T}}\in\mathbb{R}^{n_{w}\times 1}$, is the vector of well rates, $\textbf{f}=\left[(\textbf{f\,\textsuperscript{\it{p}}})^{\textnormal{T}},(\textbf{f\,\textsuperscript{\it{i}}})^{\textnormal{T}}\right]^{\textnormal{T}}\in\mathbb{R}^{n_{w}\times 1}$ and $\sum\limits_{k=1}^{n_{p}}{f}^{p}_{k}=\sum\limits_{k=1}^{n_{i}}{f}^{i}_{k}=1$ (with all elements of f between 0 and 1), where superscripts $p$ and $i$ indicate production and injection wells.

With this parameterization, the well control optimization problem can be expressed as:

[TABLE]

where J is the objective function of interest (negative NPV in our case) and $F\in\mathbb{R}^{n_{w}\times 1}\subset[0,1]$.

We address the optimization problem stated in Eq. 8 using a two-step optimization approach. This choice allows us to break the well control optimization problem into two subproblems involving different objective functions. In the first step, we optimize f:

[TABLE]

where $J_{a}$ is the first objective function. In the second step we optimize Q using $\textbf{f}^{\ast}$ from the first step:

[TABLE]

where $J_{b}$ is the second objective function.

The objective function used in the determination of optimal well-rate ratios corresponds to a measure of velocity variability. By minimizing this quantity, a more uniform sweep is achieved. This minimization is formulated as a QP problem, under the assumptions of unit mobility ratio ($M=1$) and VRR = 1. Given the optimal well-rate ratios $\textbf{f}^{\ast}$, $Q$ is determined in the second subproblem, by running a limited number of full-physics simulations, such that the NPV of the actual problem is maximized. We now describe the two subproblems in turn.

3.1.1 First subproblem: optimization of f

We introduce significant approximations, in the first step of the ST, relative to the (target) optimization of the full-physics problem. Specifically, our goal in the first step is to minimize a measure of velocity variation within the model, which we expect to lead to improved sweep and thus higher oil production and NPV (this expectation will be confirmed below). In this step we assume unit mobility ratio, which results in a linear flow equation. We note that other simulation-based approaches such as flow diagnostics [31, 30] could also be used in this step, though it is not clear that such approaches would provide problems amenable to global optimization.

The $M=1$ pressure equation can be recovered from the two-phase flow system (Eqs. 1 and 2) by setting $\mu_{w}=\mu_{o}=\mu$, and the relative permeabilities equal to their corresponding phase saturations; i.e., $f_{w}(S_{w})=k_{rw}=S_{w}$ and $f_{o}(S_{o})=k_{ro}=S_{o}=1-S_{w}$. These assumptions simplify the governing equation, which in the case of slightly compressible systems can be expressed as:

[TABLE]

where $c$ is compressibility. Darcy velocity in this case is given by:

[TABLE]

In our implementation, we solve Eq. 11 numerically using a standard finite volume implementation. Velocity components, constructed from this pressure solution, are used in the first-stage optimization. For a 3D Cartesian system, with blocks indexed as $i$, $j$ and $k$ in the $x$, $y$ and $z$-directions, and corresponding block interfaces designated $i\pm 1/2$, $j\pm 1/2$ and $k\pm 1/2$, the velocity components in the $x$-direction, across block faces $i\pm 1/2$, are given by:

[TABLE]

Here $A$ is the area of the interface and $(\tilde{q})_{i\pm\frac{1}{2},j,k}$ are the block-to-block flow rates. These are given by, e.g., $(\tilde{q})_{i-\frac{1}{2},j,k}=(T_{x})_{i-\frac{1}{2},j,k}\Delta p$, where transmissibility $(T_{x})_{i-\frac{1}{2},j,k}=\frac{k_{h}A}{\mu\Delta x}$, where $k_{h}$ is the weighted harmonic mean of $k_{x}$ in blocks $i,j,k$ and $i-1,j,k$, $\Delta x$ is the center-to-center distance between these two blocks, and $\Delta p$ is the pressure difference between the two blocks. Analogous expressions provide the fluxes and velocity components across the other grid block faces.

The velocity components at grid block centers can be obtained by averaging the velocity components at the faces:

[TABLE]

We denote the velocity components for grid block $i,j,k$ as $(u_{d})_{i,j,k}$, with $d=x,y,z$.

The linearity of Eq. 11 implies that both the pressure and the velocity fields can be expressed as a superposition of the responses of individual wells [41, 42]. This also allows the use of response matrices for pressure and velocity computations [15, 41], which can accelerate considerably the calculation of such quantities when well rates are changed. Specifically, the velocity response matrices $\textbf{V}_{d}\in\mathbb{R}^{n_{b}\times n_{w}}$, where $n_{b}$ is the number of grid blocks, enable the resulting velocity fields $\textbf{\^{u}}_{d}\in\mathbb{R}^{n_{b}\times 1}$ to be written as linear combinations of the specified well rates q̂:

[TABLE]

where $d=x,y,z$.

The methodology to build the response matrices consists of running a number of single-phase flow simulations to obtain a series of velocity coefficients. These coefficients are simply the velocity field components at every grid block of the model resulting from a single well producing at a unit flow rate. For well control optimization problems, where the well number, type, and locations are fixed, the number of single-phase flow simulations required is equal to the number of wells $n_{w}$. In each simulation, one of the wells operates individually, at a constant flow rate, until the reservoir pressure reaches pseudo-steady state (PSS); i.e., until $\frac{\partial p}{\partial t}$ becomes a constant. When this occurs, the velocity components are saved as columns in the velocity response matrices.

It is also possible (and faster) to compute the PSS solution directly rather than through the time-integration of Eq. 11. The direct PSS solution can be computed by replacing the right-hand side of Eq. 11 by $\phi K$, where $K$ is an arbitrary constant. We use time integration here since AD-GPRS is applied for these computations, and this simulator does not currently provide the PSS solution directly. We note finally that there are other ways to construct the velocity response matrices. For example, since injection and production balance in our problem, we could use $n_{w}-1$ runs, each involving two of the wells (instead of the $n_{w}$ single-well runs used in this work). The construction of the response matrices is fast, however, and the additional PSS solution required by our approach has very little impact on total elapsed time.

Figure 1 depicts the assembly of the velocity response matrix for a $3\times 2$ model with wells located at grid blocks 2 and 4. The two model setups in Fig. 1a,b indicate the well configurations for the two simulation runs required to build the response matrices. We denote $(v_{d})_{ij}$, an element of $\textbf{V}_{d}$, to be the velocity component for grid block $i$ due to the operation of a single well located in grid block $j$. The process is analogous for 3D cases.

Even though the $\textbf{V}_{d}$ matrices are constructed from pseudo-steady state solutions considering one well at a time, the multiwell response obtained through superposition will reproduce the steady-state solution as long as $\textnormal{VRR}=1$. This is now demonstrated by comparing the velocity fields generated through direct simulation (i.e., by solution of Eq. 11 until steady state is reached) and by superposition for a $100\times 100$ model containing four production and two injection wells. The channelized model and well locations are shown in Fig. 2. The resulting velocity fields in the $x$ and $y$-directions, as well as cross plots comparing the two solutions, are presented in Fig. 3. The velocity fields constructed via direct simulation and by superposition agree closely. The average absolute errors for velocities computed through superposition are about $7\times 10^{-5}$ for both components. Although the concept and use of velocity response matrices is not new [41], we are not aware of this approach being used before to construct steady-state solutions from the superposition of individual well responses obtained at pseudo-steady state.

The fact that the velocity field can be expressed as a linear combination of the well rates, using Eq. 15, means that the minimization of the (weighted) squared velocities can be formulated as a QP problem. As indicated earlier, through use of such a minimization, we can find the set of well rates f that result in a more uniform sweep and thus (expected) higher oil recovery and/or NPV. The general QP minimization can be stated as follows:

[TABLE]

where $\textbf{W}_{d}\in\mathbb{R}^{n_{b}\times n_{b}}$ with $d=x,y,z$ is a diagonal matrix whose elements are positive weighting factors. The field rate $Q$ is arbitrary at this point, and here we set $Q=1$. In the results below we specify $\textbf{W}_{d}=\textbf{I}$ for simplicity. However, by setting the elements of $\textbf{W}_{d}$ more generally, we can emphasize velocity uniformity in particular areas of interest in the reservoir. This could be beneficial, for example, when new wells are introduced into a reservoir that has already been under production for some time.

Solution of the optimization problem in Eq. 16 provides well-rate ratios that minimize the sum of the resulting squared velocities over the entire model domain. As we will see in Sect. 4.1.1, this quantity does indeed correlate with NPV for the example considered. Furthermore, the QP problem from Eq. 16 involves a quadratic objective function subject to linear constraints, and it can be solved efficiently using well-established optimization techniques that guarantee that a global optimum is found (as discussed in Sect. 2.2). The use of the response matrix approach enables the very fast computation of ${\textbf{\^{u}}}_{d}$ (i.e., we do not need to solve any flow equations), thus accelerating the determination of $\textbf{f}^{\ast}$. Note that the minimization of the actual velocity variance, rather than the squared velocity, also corresponds to a QP problem. We implemented both of these minimizations, and found that the minimization of squared velocity provided better results for the cases evaluated. This may be because the use of squared velocity penalizes high velocities more strongly. Because high extremes in the velocity field can lead to early breakthrough and inefficient displacement, and thus suboptimal NPVs, the avoidance of these effects may be very beneficial.

The computational cost for constructing the response matrix and solving the optimization problem in Eq. 16 is around 30–40% of a single full-physics AD-GPRS simulation for the cases considered here. However, as noted earlier, the PSS problem could be solved directly, without integrating in time. Such a treatment, though not implemented in this work, will decrease the cost of this step to just a few percent of a full-physics AD-GPRS run.

We reiterate that the optimization problem posed in this first ST step is heuristic. The problem does not include many of the quantities that contribute directly to NPV, such as oil price, water costs and discount rate. In addition, this problem entails simplified flow physics. Thus it is useful to assess the correspondence between minimizing squared velocity (Eq. 16) and maximizing NPV. This will be considered in Sect. 4.1.1.

3.1.2 Second subproblem: optimization of $Q$

Once the optimal well-rate ratios $\textbf{f}^{\ast}$ have been determined, the actual well rates become a linear function of the field rate $Q$, as indicated by Eq. 7. Thus, the optimization now involves only a single variable. At this stage we consider the actual objective function (in this case NPV as defined in Eq. 5) and we simulate the full-physics system.

The second-stage optimization problem can be written as:

[TABLE]

At this stage we do not require $\textnormal{VRR}=1$ or the assumption of slight compressibility, though it is expected that the overall ST will perform better when these conditions are approximated. Since the well-rate ratios $\textbf{f}^{\ast}$ are specified, $\textnormal{VRR}=1$ will still be naturally enforced unless other constraints, such as well BHP limits (discussed below), become active. We treat the problem defined in Eq. 17 very simply. Specifically, we define a range for $Q$ such that, over the simulation time frame, between 0.5 pore volumes injected (PVI) and 2.5 PVI are introduced into the model from all of the wells. We compute even increments in PVI based on $n_{\textrm{proc}}$, where $n_{\textrm{proc}}$ is the number of available processors (or the number of available simulator licenses if we have fewer licenses than processors), and then simulate each $Q$ value on a different processor. For example, for $n_{\textrm{proc}}=41$, we simulate cases corresponding to $0.5,0.55,0.60,\ldots,2.45,2.5$ PVI. This provides the green points shown in Fig. 4 in an elapsed (wall-clock) time corresponding to the time required for one full-physics simulation. We then take the $Q$ that provides the best NPV to be $Q^{*}$.

We could of course iterate this procedure further, by running additional cases in the most promising PVI regions, to provide the global optimum (red point in Fig. 4). Given the level of approximation inherent in the overall surrogate treatment, however, this additional iteration may not be necessary, assuming a sufficient number of processors (and simulator licenses) are available.

We note finally that in the special case where the system actually corresponds to $M=1$ (or is very close to $M=1$), it is possible to use highly efficient approaches such as the time-of-flight computations implemented in the Matlab Reservoir Simulation Toolbox, MRST  [43] in place of full-physics simulations. In such cases it will be very inexpensive to perform an additional iteration in the second stage of the ST.

At this point it is useful to review the treatments and approximations introduced into our two-stage ST. In the first stage, we minimize a measure of velocity variability under the assumption of $M=1$ and $\textnormal{VRR}=1$. The objective function used in this minimization is not the actual objective function, though it does correspond to a quantity that is expected to impact displacement efficiency. Then, given the optimized well-rate ratios from this first step, in the second step we maximize the actual objective function (NPV) using the full-physics model. Global optima can be determined for both problems – the first-stage optimization is a QP problem, and the second-stage optimization involves only a single variable. Since the second-stage problem can be solved in an elapsed time corresponding to a single full-physics run, and given that the first-stage problem can be solved even faster, the elapsed time for an optimization using our ST is the time required for 1–2 full-physics runs.

Another aspect of our two-step ST is that there is not any iteration between the first and second steps; i.e., each problem is solved only once. In many two-step optimization procedures, there is iteration between the two steps. In such cases, convergence to a local minimum is often assured. The lack of iteration renders our approach very efficient, but we cannot claim that it convergences to a local minimum of the actual problem. In Sect. 3.2, we will describe a treatment for BHP constraint handling that does entail iteration between the two steps, but this procedure would only be used if a BHP constraint violation is detected. We note finally that the settings provided by the two-step ST can be used as an initial guess for a traditional optimizer (e.g., MADS), and in this optimization we do expect to converge to a local minimum of the actual problem.

3.2 BHP constraint handling techniques

The two-step ST described above did not include the treatment of any constraints. In the case of well-rate specifications, as are applied here, the most important constraint is expected to be a BHP constraint. We now describe two heuristic approaches for handling this constraint. The first approach (which is the one actually implemented in this study) is applicable when the constraint violations involve a relatively small fraction of wells. The second approach, which is more complicated and involves iteration between the two ST steps, might be required if a large fraction of the wells are under BHP control.

We are concerned mainly with cases for which $\mu_{o}\gtrsim\mu_{w}$, and thus $M\gtrsim 1$. In such cases, less viscous (more mobile) fluid is injected into the reservoir, so the differences in pressure between injectors and producers tend to decrease with time. Therefore, if BHP constraint violations do occur, this will often be at the start of the simulation (we check for BHP constraint violations just after the initial transient period, which lasts for a few time steps). These constraint violations will be detected in the second stage of the ST, where full-physics simulations are performed to determine $Q^{*}$. In such cases, we specify the BHP for the affected well(s), at the maximum or minimum, and then reapportion the remaining injection and production. Specifically, well-rate ratios in this case are assigned as:

[TABLE]

where $(f^{\prime})^{p}_{k}$ is the updated well-rate ratio for production well $k$, and $n_{pr}$ is the number of production wells that operate under rate control. An analogous procedure is applied to the injection wells when a BHP constraint violation occurs at the start of the simulation.

The well-rate distributions now correspond to rescaled $\textbf{f}^{\ast}$ for those wells that remain under rate control. BHP constraints are satisfied for the other wells. We again vary rate $Q$ and compute NPV in the second ST stage. The resulting curve, which differs from that of the unconstrained case, is shown in Fig. 5. Note that the optimal rate $Q^{\ast}$ calculated in this step is not the actual total rate of the field. Rather, it is the total rate of the wells operating under rate control. The total rate of the field is the sum of $Q^{\ast}$ and the rate injected or produced by the wells operating under BHP control. On the other hand, this treatment does not impact the efficiency of the second-stage optimization – it can still be achieved, in an elapsed time of 1–2 full-order simulations, as described above.

When Eq. 18 is applied, the total rate of the field will not, in general, be strictly constant as a function of time. Field injection and production rates do, however, quickly ‘equilibrate,’ meaning VRR=1 is again achieved (and sustained throughout the simulation). This typically occurs because another well also reaches a BHP constraint, and the specification of BHP for two wells, and rates for the other wells, results in balanced injection and production in a nearly incompressible system. Nonetheless, the $\textbf{f}^{\ast}$ values determined in the first ST step for the wells under BHP control are no longer honored. This acts to reduce the benefit afforded by minimizing the squared velocity in the first ST step. However, in cases where the fraction of BHP-controlled wells is relatively small, the overall impact on the two-step ST result is not expected to be large. We reiterate that the optimal rate determined in the second step fully accounts for wells operating under BHP control.

BHP violations may also occur later in the simulation run (i.e., after the early transient period, which lasts for a few time steps), though this is less of an issue for $M\gtrsim 1$, as noted earlier. If this happens, the well is switched to BHP control, and it remains on BHP control for the rest of the run. No additional corrective measures are taken, and the other wells are kept on rate control, unless they reach their BHP limit.

Wells may also be specified to inject or produce under BHP control for operational reasons. In addition, under cases with high discount rate, early production will be preferred so rates may be higher. This could potentially lead to a larger fraction of wells reaching their BHP limits. Thus a more formal approach for handling BHP constraints may be required. We now describe such a treatment. Note that we have not implemented this approach since BHP constraints are not often reached in our runs, and the heuristic treatment described above suffices for present purposes.

For cases where a large fraction of wells reach their BHP limits, or are specified to be under BHP control for operational reasons, we believe the following iterative procedure will be effective. In the first iteration, the two-step ST is performed exactly as described earlier. In the second ST step, however, after the early transient, we compute the maximum values of $f^{i}_{k}$ or $f^{p}_{k}$ corresponding to wells that reach their BHP limits. Then, the first ST step is repeated, with the BHP-constrained $f^{i}_{k}$ or $f^{p}_{k}$ specified as upper bounds, for the corresponding wells, in the QP optimization. This will ensure that the rate allocation is indeed optimized (under the reduced-physics assumptions in the first ST step) for the BHP-constrained case. The second ST step is then repeated to determine $Q^{\ast}$ given the new $\textbf{f}^{\ast}$ values. This entire procedure can be repeated (iterated) as necessary to assure no wells violate BHP constraints for the final $\textbf{f}^{\ast}$ and $Q^{\ast}$. We note that this iterative treatment directly links the two ST steps, since the results from second step impact the first step.

A closely related issue is the treatment of wells that, for practical reasons, do not operate under rate control for extended periods, as discussed by Ushmaev et al. [14] . The strategy for handling such wells would be (essentially) that described above for wells that reach their BHP limit. In the case of wells specified to operate under BHP control, we would again iterate through the two-step ST, with upper-bound constraints for $f^{i}_{k}$ or $f^{p}_{k}$ determined from the second ST step. In this way the resulting $\textbf{f}^{\ast}$ and $Q^{\ast}$ would account for the fact that one or more wells are specified to operate under BHP control.

Practical problems may also involve nonlinear output constraints. Examples of these include the specification of maximum well water-oil ratio, or minimum well oil production rate. In general, such constraints can be handled both within the optimizer (using, e.g., a filter method as described in [38]), or directly in the simulator using well-control rules. We do not consider these types of constraints in this work. However, since the second ST step involves full-physics simulations, it would be straightforward to include these nonlinear constraints in the simulator. It might even be possible to account for them (approximately) in the iterative procedure described above.

3.3 General workflow

The two-step ST formulation entails the assembly of the velocity response matrices followed by the two optimization procedures. The Stanford Smart Fields Unified Optimization Framework (UOF) [44, 45] is used to test and compare the different optimization methods and treatments.

As explained earlier, for the construction of the velocity response matrices, a single-phase flow simulation is performed for each well, with the well specified to operate at constant flow rate. These simulations are run until the reservoir reaches pseudo-steady state. Once this is achieved, the pressure field is recorded, Darcy velocities are computed, and the velocity response matrices are assembled.

Algorithm 1 below describes the determination of ${\textbf{f}}^{*}$ and $Q^{*}$ using the two-stage ST. We first solve the squared-velocity-minimization problem (Eq. 16) under the assumption of $M=1$ and $\textnormal{VRR}=1$. Then, NPV for the full-physics problem is maximized (Eq. 17). The QP problem in the first step is solved in Matlab [46] using CVX, a package for solving convex optimization problems [47]. The second optimization step requires running $n_{\textrm{proc}}$ full-physics simulations, where $n_{\textrm{proc}}$ is the number of computing nodes or simulator licenses available. If the full-physics model corresponds to $M=1$, MRST can be used to accelerate the second-stage computations.

4 Example cases

In this section, we apply the two-step ST described in Sect. 3 to optimize well controls for 2D and 3D reservoir models. Performance is assessed for cases involving both single and multiple control periods. In the latter case, the single control period solution is used as the initial guess. The examples consider $M=1$, 3 and 5, where $M=\mu_{o}/\mu_{w}$. For $M$ = 3 and 5, we use the relative permeability curves shown in Fig. 6. For $M=1$, we specify $k_{rw}=S_{w}$, $k_{ro}=S_{o}$ and $\mu_{w}=\mu_{o}=1$. Thus the pressure and velocity fields for $M=1$ cases correspond exactly to those for single-phase flow. BHP constraints are imposed in these examples, but they are not violated at the start of the simulation in the cases presented. As noted earlier, Stanford’s Automatic Differentiation-based General Purpose Research Simulator (AD-GPRS) [36] is used for all simulations. Because the two-step ST is noninvasive with respect to the flow simulator, it can be used with any simulator.

Conventional optimization, using adjoint-gradients with SNOPT [48] and pattern search with MADS [10], is also performed in this study to enable comparison with ST results. These procedures are only guaranteed to find local optima, with the result depending on the initial guess, though in many cases the different local optima that are found correspond to very similar objective function values. This is not always the case, however. For example, in a recent study involving compositional models with nonlinear constraints, Kim [44] found local optima that corresponded to very different objective function values. In this study we consider oil-water problems subject to only bound constraints, though we do observe variation in the objective function values, in some cases, using conventional methods. For this reason, five runs with different initial guesses are performed for the adjoint-gradient and MADS optimizations whenever these methods are applied.

In this study, adjoint-gradient optimizations are terminated when the gradient satisfies a threshold of 0.001, and the optimization converges to a local optimum, or after a maximum of 15 major iterations have been performed. MADS is terminated if the stencil size decreases beyond 0.01 (where this value corresponds to the scaled controls), or after a maximum of 15 MADS iterations (for 2D models) or 20 MADS iterations (for 3D models) have been performed. Each MADS iteration requires $2n_{\textrm{opt}}$ function evaluations (full-physics simulations in our case), where $n_{\textrm{opt}}$ is the number of optimization variables [10]. Alternative stopping criteria for MADS, such as those introduced by Boch [49], could also be considered.

4.1 2D reservoir model

Our first example involves the 2D reservoir model shown in Fig. 2 (this model was generated by [37]). This model, of dimensions $100\times 100$, represents a complex channel-levee shale fluvial system. Relevant simulation and optimization parameters are provided in Table 1. In our first set of examples, only a single control period (with well rates as the control variables) is used, so we have $n_{\textrm{opt}}=n_{w}-1$.

4.1.1 Results for different well scenarios

Although field development optimization is not the main focus of this work, we first illustrate the potential use of the two-step ST within this context. This enables us to assess ST performance for a wide range of well configurations. These configurations are those actually generated during the course of a field development optimization run. Specifically, we ran the particle swarm optimization (PSO) algorithm described in [37], with a PSO swarm size of 40. The goal was to maximize NPV by optimizing the number, type, locations and controls of a set of wells, with $M=3$. We extracted the configurations (well locations and types) corresponding to the 40 particles at a relatively early stage of the optimization – specifically at iteration 25 (an early iteration was considered to maintain diversity within the swarm). We then applied ST for each of the 40 configurations. ST results are compared to optimization using MADS for the same sets of wells.

Figure 7 displays a cross-plot of the NPVs obtained by ST and MADS. Each point represents one of the 40 PSO well scenarios. Relatively close correspondence between the ST and MADS NPVs is observed (correlation coefficient $R^{2}=0.94$), with most points falling near the 45-degree line. These results demonstrate that the two-step ST is able to generate reasonable optimization results, over a range of NPVs, for the particular specifications in Table 1. In other words, the two-step ST is able to perform satisfactorily for either poor, average or favorable sets of well types and locations. This is an important capability if the method is to be used for field development optimization, in which case a large fraction of the proposed well locations are far from optimal, and it is important that they be assessed accurately.

We now perform further assessments for some of the well configurations considered in Fig. 7. These scenarios are represented by the red points in Fig. 7, and the corresponding well configurations are shown in Fig. 8a,b,c. It is evident from Fig. 7 that these cases cover a broad NPV range; i.e., we do not consider only the most promising well configurations.

In the first assessment, well scenario 3 (Fig. 8c) is considered. Our intent here is to evaluate the correlation between NPV and the sum of the squared velocities over the entire domain (this is the quantity minimized in the first ST step). In this evaluation, we specify $M=1$ and VRR = 1. In addition to the optimal well-rate ratios computed in the first ST step (${\textbf{f}}^{*}$), we consider 39 additional f vectors. One of these vectors corresponds to equally distributed rates, and the other 38 are generated randomly. In all cases the $f_{k}^{p}$ and $f_{k}^{i}$ sum to unity. For each of the 40 f vectors, the optimal field rate (${Q}^{*}$) is computed using the second ST step. This ${Q}^{*}$ determination is performed (separately) for discount rates $d$ of both 0 and 0.1 (10%).

Figure 9 displays plots of NPV versus the corresponding sum of the squared velocities for each of the 40 f vectors, for $d=0$ and 0.1. The dark blue point in each plot corresponds to the two-step ST result, the red point to the equally distributed rates case, and the light blue points to the 38 cases with random f vectors. The scale for the $x$-axis is arbitrary, since the results of the first ST step do not depend on $Q$ (here we use $Q=1$). In both plots we see that NPV is indeed inversely correlated with the sum of the squared velocities, and that the maximum NPV corresponds to ${\textbf{f}}^{*}$. These results support the choice of the objective function used in the first ST optimization step (Eq. 16). It thus appears that, even though many physical aspects of the problem are not included in the first ST step, the method can nonetheless provide useful results for well-rate ratios, at least for this example.

It is also noteworthy that the ${\textbf{f}}^{*}$ solution provides the maximum NPV with $d=0.1$ (Fig. 9b). This suggests that minimizing the squared velocities is a reasonable treatment even with some amount of discounting. We believe there are cases, however, for which this direct correspondence will no longer hold. For example, in a problem with large $d$ (e.g., $d\gtrsim 0.2$), we would expect early production to be very highly valued, and this effect might not be captured in the first ST step. Other economic aspects, such as very high injected water cost, might also represent important effects not captured in the first ST step. It is possible, however, that some of these effects could be (at least partially) treated by modifying the first-step optimization problem, or by adjusting the ${\textbf{W}}_{d}$ matrices in Eq. 16.

Next, in order to assess the degree of improvement achieved by the two-step ST (and by MADS and adjoint-gradient-based optimization), we evaluate NPV for three base cases for each of the well scenarios shown in Fig. 8. The base cases entail equally distributed injection and production rates; i.e., $f_{k}^{p}=1/n_{p}$ for $k=1,\ldots,n_{p}$ and $f_{k}^{i}=1/n_{i}$ for $k=1,\ldots,n_{i}$, for field rates that cover a range of PVI values. We ensure that the range includes the $Q^{*}$ value provided by ST.

The base-case NPVs are presented in Table 2. The $Q^{*}$ value for each case is also given. Although the use of equally distributed rates may seem somewhat simplistic, we did confirm that, for the best base cases, the water-oil-ratio (WOR) in the production wells does not reach the point where water production cost exceeds oil revenue. This means that the use of a standard reactive control procedure (in which production wells are shut in if they become uneconomic) in these cases would not have any impact on NPV. Thus, the best base cases are reasonable in this regard.

The base case results can then be compared to those in Table 3, where we show optimization results for each configuration for three different mobility ratios. This table includes ST results along with those using MADS and adjoint gradients. Results for the latter two procedures are median NPVs from five separate runs. The minimum and maximum NPVs for these methods are shown in Table 4. We reiterate that ST provides (essentially) global optima in both steps, so it does not need to be run multiple times to avoid poor local optima.

Base case and ST results are also compared, in terms of NPV versus time plots for each well configuration (for $M=3$), in Fig. 8d,e,f. These results, along with those in Tables 2 and 3, demonstrate the clear improvement provided by our two-step ST. This improvement is due both to the optimization of well-rate ratios and to the optimization of $Q$, with the relative impact of these two steps varying from case to case. Of essential interest is the close correspondence between ST optimization results and those from MADS and adjoint gradients (Table 3). The largest relative discrepancy is for well scenario 2, $M=5$, where we see ST underperform MADS by 5.6%. The close correspondence between results from our two-step ST and those from formal optimization methods is very encouraging, and suggests that the surrogate treatments may indeed capture important aspects of the true full-physics optimization problem.

The results in Table 4 demonstrate that MADS achieves very similar results in all five runs for each case. This is also generally observed for the adjoint-gradient procedure, though the well scenario 3, $M=5$ case is an exception. For the minimum NPV (1.11 BN $) run, we applied a tighter stopping criterion and performed additional SNOPT iterations, but the NPV did not improve. Thus, this result may correspond to convergence to a relatively poor local optimum.

Elapsed (wall-clock) times for the optimizations are shown in Table 5. These values are median timings for the five MADS and adjoint-gradient runs (recall ST need only be run once). Individual full-physics simulation runs for these cases require $\sim$40–50 seconds. We see that ST, as noted earlier, entails an elapsed time corresponding to 1–2 full-physics runs. MADS requires five or more iterations for these cases (we verified that additional iteration with tighter termination criteria did not improve the MADS results), so the speedup using ST is fairly significant (speedup factor of $6\times$ or more). Speedup relative to adjoint-gradient-based optimization is larger because this procedure does not readily parallelize, in contrast to MADS.

The results in this section demonstrate that our two-step surrogate treatment provides reasonable accuracy relative to formal optimization methods, and that it can lead to significant speedup. As such, it may be useful for both standalone well control optimization applications, and as an inner-loop in field development optimization. The latter point, highlighted by the ST results for many different configurations (Fig. 7), will be exploited in a future study. We next consider the application of ST for a 3D example.

4.2 3D reservoir model

The 3D model used in this assessment is shown in Fig. 10. This $40\times 30\times 9$ model (10,800 total grid blocks) represents a prograding fluvial channel system, and corresponds to a portion of a synthetic geological model developed by Castro [50]. We assess ST performance for two different scenarios involving (1) three producers and three injectors, and (2) four producers and two injectors. All wells are completed over the full reservoir thickness. The two scenarios are shown in Fig. 11. Both single and multiple control periods are considered for this case.

Relevant simulation and optimization parameters appear in Table 6. We again consider $M=1$, 3 and 5, and use the same relative permeability curves as in the 2D cases (the curves in Fig. 6 for $M\neq 1$, and straight-line relative permeabilities for $M=1$).

4.2.1 Single control period

Optimization results for the two scenarios shown in Fig. 11 will now be presented. We first consider the case where well rates are fixed for the entire simulation time frame (as was the treatment in all of the 2D runs). As in the 2D case, several base cases, with equal distribution of rates among all injection and all production wells, are assessed for $M=3$. Again, for the best base cases, none of the production wells reach the point where water production cost exceeds oil revenue. As noted earlier, this indicates that the use of standard reactive controls would not impact the best base-case NPVs.

The base-case NPVs are shown in Table 7, and the optimized NPVs using ST and MADS are presented in Table 8. The median, minimum and maximum NPVs using MADS are provided. Adjoint-gradient-based optimizations were not performed for this case, though we would expect these results to closely coincide with the MADS results, as in the 2D example.

As in the 2D cases, we again see that the NPVs provided by ST clearly exceed those for the three base cases for both scenarios. Base-case and ST results for cumulative field oil and water production and cumulative water injection, for well scenario 2 with $M=3$, are shown in Fig. 12. Results are also presented for total injection (and production) of $Q^{*}$, as determined by ST, but with equal injection and production between wells. The ST results display the most oil production, along with the least water production, of the cases considered. These results highlight the importance of optimizing both f and $Q$.

Comparison of ST and MADS results in Table 8 demonstrates that ST again provides results quite close to those achieved using a formal optimization method. The largest relative difference between ST and (median) MADS results is for well scenario 1, $M=5$, where ST underperforms MADS by 4.7%. For well scenario 1, $M=1$, ST actually achieves a larger NPV than any of the MADS runs. We also see that there is very little variation between the MADS results in most cases. For well scenario 1 with $M=1$ and 3, the MADS results do vary by 1–2% between the minimum and maximum NPV runs. For the minimum-NPV runs, we tightened the termination criterion and performed additional MADS iterations, but this did not provide improvement in the NPVs. This observation supports the idea that the different local optima may correspond to slightly different NPVs.

Table 9 presents elapsed time, under full parallelization, for ST and MADS optimizations. Individual simulation runs require $\sim$105–115 seconds for these cases. We thus see that ST elapsed times continue to be less than the time required for two full-physics runs. MADS requires six iterations for these cases, and speedup using ST of about $5\times$ relative to MADS is observed. It is again apparent that ST is able to achieve reasonable results along with substantial speedup.

4.2.2 Multiple control periods

Our examples up to this point have involved a single control period and zero discount rate (except for the $M=1$ results in Fig. 9b, where we considered $d=0.1$). We now relax these assumptions and consider multiple control period strategies, with and without discounting. Both cases involve well scenario 1 (three producers and three injectors). In the first case we consider $M=3$ and $d=0$, and in the second case $M=5$ and $d=0.05$. The well rates are now determined four times during the simulation time frame, which increases $n_{\textrm{opt}}$ to 20 (16 well-rate ratios plus 4 field rates). Two approaches can be used to assess ST performance relative to MADS for this case. First, we can simply apply the ST well rates over the entire time frame. Second, the single-period ST solution can be used as an initial guess for MADS (we designate this case ST & MADS).

The NPV results and corresponding elapsed times obtained with ST (standard single control period solution), ST & MADS, and MADS (we give the median, minimum and maximum NPVs for the five standalone MADS runs) are presented in Tables 10 and 11. We see that both methods provide generally similar NPV results. For the first case, there does not appear to be substantial benefit from using multiple control periods (compare results in Tables 8 and 10). We do, however, observe that the NPV resulting from the ST & MADS approach provides an objective function value that is 3% higher than the median solution from standalone MADS for the first case. Even though this difference is relatively small, this result suggests that the ST solution is indeed a reasonable initial guess for the subsequent MADS run. The corresponding optimized injection and production rates for the four control intervals, obtained by ST & MADS for the $M=3$, $d=0$ case, are presented in Figs. 13 and 14. We do see some shifts from control period to control period, though these shifts are relatively slight.

For the second case ($M=5$, $d=0.05$), the ST result is 2.8% below the median MADS result, but it is above the minimum MADS result. The ST & MADS approach provides the same objective function value as the median MADS case. These results, along with those in Fig. 9b, suggest that our two-step ST is indeed applicable for cases with nonzero discount rate. As noted earlier, the first ST step may become less reliable in cases with large $d$, where very early production is highly beneficial. The effect of large $d$ would, however, be captured in the second ST step.

From the timings in Table 11, we see that ST & MADS is slightly faster than MADS (ST & MADS requires 70–80% of the median time for MADS). The improvement in NPV with iteration for the first case ($M=3$, $d=0$) is shown in Fig. 15, which displays the evolution of NPV for ST, ST & MADS, and standalone MADS (median case). In this case, ST & MADS leads to speedup relative to standalone MADS because it requires 9 iterations, while standalone MADS performs 13  iterations (both methods use the same termination criteria).

A potential advantage of ST & MADS is that standalone MADS might be run multiple times to find the best solution (even though the differences in NPVs between MADS runs are typically rather small), while ST & MADS needs to be run only once. This is possible because the initial guess provided by standard ST is expected to lead to a promising (local) optimum, as is indeed observed for this example. This behavior, if achieved consistently over a range of more challenging test cases, would represent a real advantage for ST & MADS, as multiple optimization runs would not be required.

5 Concluding remarks

In this work, we introduced a new two-step surrogate treatment (ST) for well control optimization of oil fields under waterflood. The method provides injection and production rates for all wells for a single control period. Our ST proceeds by first determining appropriate well-rate ratios, which are expected to lead to more uniform reservoir sweep. In the second step, overall field rate is optimized such that the objective function of the actual problem (NPV in the cases considered here) is maximized or minimized. In the first step we use simplified physics, specifically the assumption of unit-mobility-ratio displacement and exact voidage replacement. This step is accelerated through use of velocity response matrices, which can be constructed efficiently in a preprocessing step. Full-physics simulations are used in the second ST step, but this minimization can be performed quickly since it involves only a single optimization variable. Near-global optima are obtained in both steps, though this does not provide the global optimum to the actual problem due to the approximations introduced in the ST. The overall two-step procedure requires an elapsed time corresponding to the runtime for 1–2 full-order simulations, assuming full parallelization in the second ST step.

The ST was applied for both 2D and 3D optimization problems. For the 2D case, 40 well scenarios (corresponding to a wide range of well locations and types), generated by PSO during a field development optimization run, were used to compare ST and MADS. Relatively close correspondence between the two sets of optimization results, for $M=3$, was observed. More detailed comparisons were presented for three of the well scenarios for systems with $M=1$, 3 and 5. Multiple MADS and adjoint-gradient runs were performed to assure avoidance of poor local optima with these methods (though this is generally not a major issue with oil-water well control optimization problems of the type considered here). Again, close correspondence between results from our two-step ST and the formal optimization procedures was observed. Similar observations applied for the 3D case, where two well scenarios were considered, again for $M=1$, 3 and 5. In terms of timings, under full parallelization ST provided speedups of $5\times$ or more relative to MADS for the cases considered. For a case involving multiple control periods, where the ST result provided the initial guess for MADS, a speedup of less than $2\times$ was observed. However, for such cases standalone MADS would likely be run multiple times, while the ST & MADS approach need be run only once.

There are a number of areas in which future work extending our two-step ST could be performed. More general well control optimization setups, involving, e.g., additional cases with nonzero discount rates or the consideration of optimal economic project life [51], could be considered. Optimization under geological uncertainty is another important aspect that should be addressed. In this case realization selection, using procedures described in [52], could be applied. The two-step ST may be particularly useful within the context of field development optimization, where well number, locations, types and controls must be determined. Such MINLP problems are extremely demanding computationally, so surrogate treatments can have a large impact. The applicability of the two-step ST for this problem is suggested by the results shown in Fig. 7, and our findings in a followup study indeed confirm this. The two-step ST could also be extended to handle more complex production scenarios, such as primary production followed by waterflood, or water-alternating-gas processes. Further testing with realistic 3D models should also be performed. Finally, comparisons with other heuristic procedures for determining well allocation/rates (such as those in [34, 14]) should be performed. Optimization results using the various approaches could then be assessed, and these findings could be used to devise a hybrid implementation based on the most promising components from the different methodologies.

Acknowledgements

We thank Nikhil Padhye for providing the Unified Optimization Framework (which facilitated testing and comparing the different optimization algorithms used in this study) and for helping us with its use. We are also grateful to Marco Antonio Cardoso for useful discussions and suggestions. We thank Petrobras and the industrial affiliates of the Stanford Smart Fields Consortium for partial funding of this work.