Dynamical systems often exhibit limit cycle oscillations (LCOs), self-sustaining oscillations of limited amplitude. LCOs can be supercritical or subcritical. The supercritical response is benign, while the subcritical response can be bi-stable and exhibit a hysteretic response. Subcritical responses can be avoided in design optimization by enforcing LCO stability. However, many high-fidelity system models are computationally expensive to evaluate. Thus, there is a need for an efficient computational approach that can model instability and handle hundreds or thousands of design variables. To address this need, we propose a simple metric to determine the LCO stability using a fitted bifurcation diagram slope. We develop an adjoint-based formula to efficiently compute the stability derivative with respect to many design variables. To evaluate the stability derivative, we only need to compute the time-spectral adjoint equation three times, regardless of the number of design variables. The proposed adjoint method is verified with finite differences, achieving a five-digit agreement between the two approaches. We consider a stability-constrained LCO parameter optimization problem using an analytic model to demonstrate that the optimizer can suppress the instability. We also consider a more realistic LCO speed and stability-constrained airfoil problem that minimizes the normalized mass and stiffness. The proposed method could be extended to optimization problems with a partial differential equation (PDE)-based model, opening the door to other applications where high-fidelity models are needed.