Simulated Annealing

Key Takeaways

• The Metropolis algorithm that underlies simulated annealing defines the acceptance probability for each proposed move in the parameter space: at each iteration, the algorithm proposes a random perturbation to the current solution (a change to one or more model parameters), evaluates the objective function at the new point (the "energy" of the system in the annealing metaphor), and accepts the new solution with probability P = 1 if the objective function improves (the new energy is lower) or with probability P = exp(-delta_E / T) if the objective function worsens (the new energy is higher by an amount delta_E), where T is the current "temperature" of the algorithm; at high temperature T, even large objective function increases are accepted with moderate probability, allowing the algorithm to explore the parameter space broadly and escape from local minima; as T decreases according to a cooling schedule (typically geometric cooling, T_k+1 = alpha x T_k with alpha = 0.90-0.99, or logarithmic cooling, T_k = T_0 / log(1 + k)), the acceptance probability for worsening moves decreases, and the algorithm focuses increasingly on local refinement around the best solution found so far; the quality of the final solution depends critically on the cooling schedule: too fast cooling replicates the failure of gradient descent methods by trapping the algorithm in local minima; too slow cooling wastes computational resources on unnecessary exploration; optimizing the cooling schedule for a specific application often requires experimentation or theoretical analysis of the objective function landscape.
• Reservoir model history matching using simulated annealing adjusts the spatial distribution of reservoir properties (permeability, porosity, fault transmissibility, aquifer strength) by iteratively modifying the reservoir model and comparing the simulated production history (pressure, GOR, WOR, and rate from each well) with the observed production history from the field; the objective function quantifies the mismatch between simulated and observed history, and simulated annealing minimizes this mismatch by exploring the space of possible property distributions; the advantage of SA over gradient-based history matching methods is its ability to explore non-unique solutions (many different property distributions may produce equally good matches to the production history) and to escape from local minima in the objective function that correspond to property distributions that match some observed data but not all; the output of SA history matching is typically not a single best-fit model but an ensemble of models that all provide acceptable matches, which can be used to quantify the range of reservoir predictions (production forecasts, injection response, drainage patterns) consistent with the observed data; the Ensemble Kalman Filter (EnKF) and Ensemble Smoother with Multiple Data Assimilation (ES-MDA) are alternative ensemble-based history matching methods that have largely replaced SA in industrial practice for large, complex models, though SA remains useful for smaller problems and for problems where the objective function is highly non-Gaussian.
• Geostatistical simulation using simulated annealing generates realizations of reservoir property fields (porosity, permeability, facies) that reproduce a target statistical distribution (histogram) and a target spatial correlation structure (variogram or training image) while honoring data values at well locations: the algorithm starts with a random initial realization and iteratively swaps values between pairs of grid cells, accepting each swap if it improves the match to the target statistics (or rejecting it with a probability governed by the SA temperature schedule if it worsens the match); this approach, called Sequential Simulated Annealing (SSA) or Simulated Annealing (SA) for geostatistics, produces geologically plausible realizations that can be used in reservoir simulation to quantify uncertainty in production forecasts; SA-based geostatistical simulation has been largely superseded in industry practice by Sequential Gaussian Simulation (SGS), Sequential Indicator Simulation (SIS), and Multiple Point Statistics (MPS) methods that are computationally more efficient for generating large 3D property realizations, but SA retains advantages for conditioning to complex or multiple data types (seismic amplitudes, well test permeabilities, production data) simultaneously.
• Well placement optimization using simulated annealing searches the space of possible well locations and trajectories to identify the development plan (number of wells, well types, and positions) that maximizes the economic value of the field over its producing life: the objective function typically combines NPV (net present value of the production forecast from the proposed well pattern) with operational constraints (minimum distance between wells, maximum well inclination, cost per well) and geological uncertainty (the range of production forecasts across the ensemble of history-matched reservoir models); SA-based well placement optimization runs the reservoir simulator many times (hundreds to thousands of evaluations) with different proposed well configurations, and uses the SA algorithm to navigate the combinatorially large space of possible well plans (the number of possible combinations of well locations, trajectories, and completions is astronomically large for fields with many potential well slots) while avoiding the premature convergence to locally optimal but globally suboptimal development plans that would result from simpler sequential well-by-well optimization approaches; the computational cost of running the reservoir simulator thousands of times limits SA-based well placement optimization to models with manageable simulation run times, or requires use of proxy models (fast surrogate response surfaces calibrated to the full-physics simulator) for preliminary screening before confirmatory full-physics evaluation.
• Seismic full-waveform inversion (FWI) and tomographic velocity model building both use SA or SA-like global optimization algorithms as alternatives to gradient-based local optimization approaches when the velocity model space has multiple local minima (cycle-skipping problem in FWI, where the simulated waveforms are more than half a wavelength away from the observed waveforms at the starting velocity model, causing gradient-based methods to converge to wrong local minima): SA can in principle escape cycle-skipping by accepting velocity models that temporarily increase the waveform misfit on the path to finding the global minimum, but the high computational cost of evaluating seismic wave propagation for each proposed velocity model makes SA impractical for full-scale 3D FWI; instead, SA-like approaches (simulated annealing with parallel tempering, basin hopping) are used for 2D velocity model building in particularly complex geological settings (salt proximity, overthrust structures) where traditional gradient-based FWI consistently fails to converge to the correct velocity model; the practical FWI workflow for salt-prone areas uses a combination of manual velocity model building, ray-based tomography, and gradient-based FWI with careful starting model construction to avoid cycle-skipping, supplemented by SA-like global optimization only for specific problem areas where local methods demonstrably fail.