Summary: This study proposes a novel approach for building surrogate models, in the form of reduced-order models(ROMs), for partial differential equation constrained optimization. A physics-constrained data-driven framework is developed to transform large-scale nonlinear optimization problems into ROM-constrained optimization problems. Unlike conventional methods, which suffer from instability of the forward sensitivity function, the proposed approach maps optimization problems to system dynamics optimization problems in Hilbert space to improve stability, reduce memory requirements, and lower computational cost. The utility of this approach is demonstrated for aerodynamic optimization of an NACA 0012 airfoil at \(Re = 1000\). A drag reduction of 9.35% is obtained at an effective angle of attack of eight degrees, with negligible impact on lift. Similarly, a drag reduction of 20% is obtained for fully separated flow at an angle of attack of \(25^\circ\). Results from these two optimization problems also reveal relationships between optimization in physical space and optimization of dynamical behaviours in Hilbert space.