Geometrical imperfections and residual stresses are two inevitable aspects in studying of graphene sheets (GSs) which despite of their vital influence on the forced nonlinear vibrational behavior they have not been considered in previous research works. Hence, this work is an attempt to fill this important gap for the first time. Utilizing Eringen's nonlocal plate theory, assuming von Karman nonlinear terms in strain-displacement relations and employing calculus of variations, we have derived partial differential equations (PDEs) of the displacement field. The function of initial curvature was assumed as the first mode-shape of the simply-supported rectangular plates with a variable mode-coefficient. Next, Galerkin decomposition method has been used to achieve the equation of motion and for nonlinear vibration analysis time multiple scales method (TMS) was applied. The nonlinear frequency response relation has been analytically achieved. Afterwards, comparing to molecular dynamic (MD) results presented in previous researches we validated natural frequencies of single layer graphene sheets (SLGSs) in different sizes. Furthermore, nonlinear frequency response curves and time histories based on TMS, were in an excellent agreement with numerical method results. Therefore, the influences of external pressure, nonlocal parameter, geometric imperfection and residual stress have been presented comprehensively. It was observed that both imperfection and stress intensify the softening behavior in the system. Interestingly, for a certain controllable magnitude of the imperfection, it has been illustrated that the generated softening term counteracts with the hardening due to geometrical nonlinearity and consequently the nonlinear system behaves as a linear one.