This work presents a fast and robust method for optimizing the stationary radial distribution of temperature, density and parallel current density in a tokamak plasma and its application to first-principle-based modeling of the ITER hybrid scenario. A new solver is implemented in the RAPTOR transport code, enabling direct evaluation of the stationary solution to which the radial plasma profiles evolve. Coupled to a neural network emulation of the quasi-linear gyrokinetic QuaLiKiz transport model (QLKNN-hyper-10D), a first-principle-based estimate of the stationary state of the core plasma can be found at unprecedented computational speed (typically a few seconds on standard hardware). The stationary state solver is then embedded in a numerical optimization scheme, allowing the optimization of tokamak plasma scenarios in only a few minutes. The proposed method is applied to investigate the performance of ITER hybrid scenarios at different values of total plasma current, plasma density and pedestal height and for different power contributions in a heating mix consisting of electron cyclotron and neutral beam heating. Optimizing the radial distribution of electron cyclotron current drive (ECCD) deposition, the q profile is tailored to maximize the fusion gain Q, by maximizing the energy confinement predicted through the first-principles-based transport model, while satisfying q > 1, avoiding sawtooth oscillations. It is found that optimal use of ECCD in ITER hybrid scenarios is to deposit power as close to the core as possible, while maintaining sufficient off-axis current drive to keep q above 1. Upper limits for the fusion gain Q are shown to be constrained either by minimum power requirements for the separatrix power flow to maintain H-mode or by minimum current drive requirements for q profile tailoring. Finally, it is shown that the ITER hybrid scenario operating window is significantly extended by an upgrade of the electron cyclotron power to 40 MW.