Limit cycle oscillations (LCO) can be supercritical or subcritical. The supercritical response is benign, while the subcritical response can be destructive. It is desirable to avoid subcritical responses by enforcing LCO stability through design optimization constraints. However, there is a need for a computational approach that can model this instability and can handle hundreds or thousands of design variables. We propose a simple metric to determine the LCO stability using a fitted bifurcation diagram slope to address this need. We develop an adjoint-based formula to compute the stability derivative efficiently with respect to many design variables. To evaluate the stability derivative, we only need to compute the time-spectral adjoint equation three times no matter the number of design variables in the optimization problem. 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 and demonstrate that the optimizer can suppress the instability. Finally, we consider a more realistic LCO speed and stability-constrained mass minimization problem for an airfoil. The proposed method can be extended to optimization problems with a PDE-based model, opening the door to applications where high-fidelity models are needed.