The optimization of a nuclear reactor's performance and stability are two important factors in its optimal control. Stability analysis of a nonlinear system with a high probability of oscillations and bifurcations poses a challenge for the explicit computation of eigenvalues at each time step. Such computations are computationally expensive and nonsmooth near bifurcation points. In this work, a novel optimal control scheme using a neural network is proposed for boiling water reactors (BWR), heavy water reactors (HWR), and pressurized water reactors (PWR). Feedforward neural networks are used to approximate the maximum real eigenvalue of the system's Jacobian as a smooth function of the system's states and bifurcation parameters. For the PWR and HWR systems, a smooth approximation of the hyperbolic tangent function, x/(1+∣x∣), is used on the network outputs to ensure smoothness in the objective function and constraints for IPOPT optimization. For the BWR systems, the network outputs are used without any further processing because solver stability is not an issue. The approach combines surrogate models within a dynamic optimization problem using Pyomo software, with a soft penalty function and a hard constraint function to prevent unstable solutions. The results show that a positive sign of the penalty parameters is effective for BWR and HWR systems, while a negative sign is effective for PWR systems. The proposed approach is computationally efficient and can be generalized to incorporate stability into optimal control for various nuclear reactor systems.