A Min-Max Gradient Search Method for Constrained Simulation Optimization

Abstract: Constrained simulation optimization (CSO) is a general framework for optimizing stochastic systems under performance constraints. It arises widely in practice where objective and constraint evaluations are available only through noisy simulation outputs. Compared with the unconstrained setting, the lack of accessible analytical gradients for simulation-based constraints makes it more challenging to develop efficient solution methods and establish non-asymptotic guarantees. To address this gap, we propose a novel single-loop algorithm, called min-max gradient search (MGS), which integrates a primal-dual framework with stochastic gradient estimators. Unlike conventional stochastic approximation methods based on gradient descent for solving simulation optimization problems, such as Zhou and Bhatnagar (2017) and Hu and Fu (2025), MGS performs alternating gradient descent and ascent on the primal and dual variables, which improves the objective while penalizing constraint violations. For the first time, we establish a finite-time convergence guarantee for single-loop CSO algorithms by showing that MGS converges to a stationary solution (a Karush-Kuhn-Tucker point under mild conditions) at a rate of $\tilde{O}(T^{-1/3})$, where $T$ is the number of iterations. Numerical experiments on a serial queuing system and a 2000-dimensional optimization problem demonstrate the superior performance and scalability of MGS.

• The paper presents a single-loop min-max gradient search algorithm that reformulates constrained simulation optimization into a primal-dual framework, ensuring convergence to KKT points.
• It uses Gaussian-smoothing finite differences with common random numbers and momentum-based updates to reduce gradient variance and scale effectively to high-dimensional problems.
• Empirical results show substantial improvements in simulation efficiency and convergence rates compared to penalty, Bayesian, and genetic algorithms on complex CSO problems.

Min-Max Gradient Search for Constrained Simulation Optimization

Problem Setting and Challenges

The paper addresses constrained simulation optimization (CSO) problems of the form:

$$\min_{x\in \mathcal{X}} \,\, h_0(x) = \mathbb{E}[h_0(x;\xi)] \quad \text{s.t.} \quad h_j(x) = \mathbb{E}[h_j(x;\xi)] \le 0, \quad \forall j \in \mathcal{I},$$

where both objective and constraint evaluations are only accessible via noisy simulation outputs due to the stochastic nature of the underlying system. This renders gradient-based optimization and feasibility checks particularly challenging, as no analytical forms are available for either objectives or constraints. The feasible set $\mathcal{X}$ is convex and compact, and the constraints further cut out a feasible region that is not known explicitly.

This setting is prevalent in operations research and applied fields where performance and quality constraints are enforced, but only estimated via Monte Carlo-type simulation (e.g., queuing systems, energy reliability under uncertainty, and inventory management with service-level targets). The simulation-only scenario rules out direct projection onto the constraint set and prevents the use of vanilla stochastic approximation methods. This creates a need for gradient estimates robust to simulation noise and algorithms that can jointly handle optimization and feasibility.

Min-Max Reformulation and Methodology

To circumvent the intractability of direct projection onto the feasible set, the problem is recast into a min-max Lagrangian framework, introducing dual variables (Lagrange multipliers) to penalize constraint violations:

$$\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \mathbb{E} \left[ h_0(x; \xi) + \sum_{j \in \mathcal{I}} y(j) h_j(x; \xi) \right],$$

with the dual variable $y$ constrained in a bounded and convex set $\mathcal{Y}$ for stability and tractable finite-sample analysis. This lifts the problem into a min-max optimization setting which, under suitable technical conditions, ensures that stationary points correspond to KKT points of the original problem.

Owing to the lack of analytical gradients, the method employs random-direction, Gaussian-smoothing finite-difference estimators for (sub-)gradients of both the primal and dual variables. Crucially, common random numbers (CRN) are used to evaluate perturbed and baseline simulations with shared randomness, substantially reducing the variance of the estimator, especially for small perturbation radii.

Two critical algorithmic enhancements are introduced:

1. Strong-Concavity Regularization: A small quadratic regularization term $-\frac{\mu}{2}\|y\|^2$ is added to the Lagrangian to induce strong concavity in $y$, stabilizing the max-ascent step in nonconvex/nonconcave settings, which is well known to mitigate cycling and promote convergence.
2. Momentum and Batch Averaging: Momentum-based updates, similar to Polyak’s heavy ball, are employed within both primal and dual steps to smooth noisy stochastic gradients. In addition, the use of batched, averaged estimators over $q$ independent random directions further reduces estimator variance.

The result is a single-loop Min-Max Gradient Search (MGS) algorithm, which alternates between primal descent (for $x$) and dual ascent (for $y$) without needing to solve inner-loop subproblems at each iteration. This design contrasts with prior nested-loop approaches, yielding substantial practical advantages in implementation and parameter tuning.

Non-Asymptotic Convergence Guarantees

A key contribution of the paper is the establishment of the first finite-time convergence guarantee for single-loop CSO algorithms addressing general nonconvex, simulation-only settings. Under standard smoothness, bounded variance, and compactness assumptions—as well as an appropriate choice of step size schedules and batch sizes—the MGS algorithm provably achieves:

• Convergence rate to stationarity: $\mathcal{X}$0 in the number of simulation iterations $\mathcal{X}$1 (Corollary 1).
• Simulation complexity: $\mathcal{X}$2 simulations to reach $\mathcal{X}$3-stationarity, improving over the previous best $\mathcal{X}$4 for zeroth-order min-max optimization.
• The bound is dimensionally robust: Unlike many random-search approaches whose rates degrade rapidly with increasing dimension, the rate here depends only weakly on $\mathcal{X}$5, allowing effective scaling to high-dimensional problems.

The analysis leverages the reduced variance afforded by CRN and momentum, optimizing the tradeoff between bias (due to finite-difference smoothing) and variance (from simulation noise). Notably, the variance constant for CRN is strictly smaller than for independent randomization, implying a constant-factor reduction in the effective replication budget required to maintain a given level of stationarity.

KKT Gap Control: The projected stationarity measure used in the min-max problem can be explicitly related to an approximate KKT gap for the original CSO problem, ensuring that convergence under the algorithm translates to both feasibility and near-optimality for the original constraints.

Empirical Results and Numerical Insights

The method is empirically validated on two representative problem classes:

• Serial Queuing System: Optimizing service rates across a network of queues to minimize cost under service-level constraints. The MGS algorithm shows rapid and stable convergence, markedly outperforming penalty-based methods, Bayesian optimization with constrained expected improvement (CEI), and genetic algorithms in terms of simulation efficiency and final solution quality.

  Figure 1: Performance comparison of MGS, penalty-based methods, Bayesian optimization, and genetic algorithms in a constrained queuing system, highlighting the superior convergence and solution quality of MGS.

• High-Dimensional CSO (2000-D): The MGS algorithm is tested on a synthetic 2000-dimensional polynomial-constrained stochastic optimization, where it attains feasibility and cost optimality reliably, while competitors either fail to scale or converge extremely slowly.

  Figure 2: Performance of MGS versus penalty-based methods in a 2000-dimensional constrained simulation optimization problem, demonstrating the scalability and robustness of the proposed algorithm.

Three key ablation experiments are also reported:

• Effect of CRN: When CRN is disabled, convergence speed and trajectory stability degrade markedly, verifying the analysis regarding variance reduction.
• Batched Gradient Estimation (Effect of $\mathcal{X}$6): Increasing replication per gradient estimator up to moderate $\mathcal{X}$7 improves convergence stability and speed, with diminishing returns as $\mathcal{X}$8 grows large.
• Momentum and Regularization: Removal of momentum diminishes, but does not catastrophically harm, convergence speed. Omitting strong-concavity regularization ($\mathcal{X}$9) induces instability only in certain nonconcave cases, but is theoretically critical to the convergence proofs.

  Figure 3: Stationarity measure and constraint violation for MGS under ablation of CRN, momentum, and strong-concavity regularization.

Strong Numerical Findings and Contrasts

The authors report superior numerical performance and scalability for MGS relative to leading alternatives on standard benchmarks. Notably:

• The convergence rate of $$\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \mathbb{E} \left[ h_0(x; \xi) + \sum_{j \in \mathcal{I}} y(j) h_j(x; \xi) \right],$$0 aligns tightly with practical performance.
• For both moderate and large batch sizes $$\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \mathbb{E} \left[ h_0(x; \xi) + \sum_{j \in \mathcal{I}} y(j) h_j(x; \xi) \right],$$1, MGS rapidly attains feasibility with minimal constraint violation, whereas penalty and metaheuristic methods converge much more slowly and either oversatisfy constraints (wasting resources) or fail to achieve feasibility.
• In very high dimensions ($$\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \mathbb{E} \left[ h_0(x; \xi) + \sum_{j \in \mathcal{I}} y(j) h_j(x; \xi) \right],$$2), Bayesian and evolutionary algorithms are computationally infeasible or ineffective, while MGS remains robust.

  Figure 4: Impact of varying batch size $$\min_{x \in \mathcal{X}} \max_{y \in \mathcal{Y}} \mathbb{E} \left[ h_0(x; \xi) + \sum_{j \in \mathcal{I}} y(j) h_j(x; \xi) \right],$$3 on solution stability and convergence speed in the serial queuing problem.

Implications and Future Directions

The min-max gradient search framework for CSO provides:

• A theoretically rigorous, practically implementable approach to simulation-based optimization under constraints, requiring only access to noisy simulation outputs.
• Substantial improvements in both finite-time sample complexity and runtime, making it feasible for large-scale, high-dimensional CSO problems common in industrial and scientific simulation settings.
• A template for integrating CRN and momentum with stochastic zeroth-order optimization, which may stimulate advances in black-box min-max problems, including adversarial ML and robust policy optimization.

Potential directions for further research include:

• Extension to multi-objective and non-convex constraint forms.
• Integration of adaptive batch size selection and learning-rate scheduling.
• Investigation of second-order or quasi-Newton enhancements.
• Application to reinforcement learning where simulation-based policy constraints are present.

Conclusion

This work introduces and analyzes the MGS algorithm, a single-loop, min-max gradient search approach tailored to constrained simulation optimization with noisy, black-box objectives and constraints (2606.07029). By leveraging stochastic primal-dual updates, Gaussian-smoothing finite-difference gradients, common random numbers, and momentum, MGS achieves provably efficient convergence with strong empirical performance in both moderate and high-dimensional settings. The proposed scheme broadens the range of tractable CSO problems and sets a new theoretical and computational benchmark for simulation optimization with functional constraints.