In the context of production optimization, we consider the general problem of finding the well controls that maximize the net present value (NPV) of life-cycle production, where the well controls are either the bottomhole pressure (BHP) or a rate (oil, gas, water, or total liquid) at each well on a set of specified control steps (time intervals), with the limitations on surface facility considered as nonlinear-state constraints [e.g., field-liquid-production rates (FLRs), field-water-production rates (FWRs), and/or field-gas-production rates]. If the reservoir simulation used for reservoir management has sufficient adjoint capability to compute gradients of the objective function and all state constraints, we show that one can develop a significantly more computationally efficient procedure by replacing the adjoint-enhanced reservoir simulator by a proxy model and optimizing the proxy. Our methodology achieves computational efficiency by generating a set of output values of the cost and constraint functions and their associated derivative values by running the reservoir simulator for a broad set of input design variables (well controls) and then using the set of input/output data to train a proxy model to replace the reservoir simulator when computing values of cost and constraint functions and their derivatives during iterations of sequential quadratic programming (SQP). The derivation of the equations for computing the proxy-based model that uses both function and gradient information is similar to that of least-squares support vector regression (LS-SVR). However, this method is referred to as gradient-enhanced support vector regression (GE-SVR) because, unlike LS-SVR, the method uses derivative information, not just function values, to train the proxy. Similar to LS-SVR, improved (higher) estimated optimal NPV values can be obtained by using iterative resampling (IR). With IR, after each proxy-based optimization, one evaluates the cost and constraint functions and their derivatives at the estimated optimal controls using reservoir-simulator output, and then adds this new input/output information to the training set to update the proxy models for predicting NPV and constraints. Using the updated proxies, one applies SQP optimization again. IR continues until the simulator and proxy evaluated at the latest estimate of the optimal well controls give the same value of NPV within a specified percentage tolerance and the constraints evaluated by reservoir simulator at the latest optimal well controls are such that the constraints are satisfied within some small specified tolerance. Our results indicate that proxy-based optimization with iterative resampling might require up to an order of magnitude less computational time than pure reservoir-simulator-based optimization. By comparing the results generated with an LS-SVR proxy with the GE-SVR results, we find that GE-SVR is roughly an order of magnitude more computationally efficient than LS-SVR but also provides a better approximation of a complex cost-function surface so that it is possible to locate multiple optima in cases where LS-SVR fails to identify the multiple optima.