Robust data-driven reduced-order models (ROMs) could enable near-optimal control for very high-dimensional nonlinear dynamical systems, with applications in active flow control such as relaminarizing turbulent flows and recovering from aerodynamic stall. With initial conditions far away from the desired steady state, accounting for nonlinear effects becomes crucial — yet makes solving the resulting Hamilton-Jacobi-Bellman (HJB) equation, which defines the value function over the continuous state space, computationally intractable due to the curse of dimensionality. Reduced-order models (ROMs) can help in a variety of ways, but existing methods often fail to capture relevant dynamics for the control problem . To overcome these challenges, we first approximate the infinite-horizon optimal control problem with a large finite horizon to obtain a time-invariant version of the value function. We furthermore employ an indirect method of trajectory optimization (which is feasible in high dimensions) to obtain state and costate data offline, along locally optimal trajectories for estimating state and gradient covariance matrices. This method builds on the Pontryagin minimum principle and other related work that establishes the costate (adjoint variables) provided as generalized gradients of the optimal value function satisfying the HJB equation . An oblique projection obtained by balancing these matrices with initial conditions sampled uniformly from an uncontrolled attractor is used to identify active directions in the state space along which the value function is most sensitive, and states have large variance . The oblique projection obtained is used to build surrogate models for both the value function and the optimal feedback control law, which is validated on the full order model. ROM-based state estimators are also built with these projections for closed-loop feedback control. We assess the quality of the resulting ROMs across linear and nonlinear PDE control problems, benchmarking against existing ROM approaches.