Exact solutions for the torsional bifurcation buckling of functionally graded (FG) multilayer graphene platelet reinforced composite (GPLRC) cylindrical shells are obtained. Five types of graphene platelets (GPLs) distributions are considered, and a slope factor is introduced to adjust the distribution profile of the GPLs. Within the framework of Donnell’s shell theory and with the aid symplectic mathematics, a set of lower-order Hamiltonian canonical equations are established and solved analytically. Consequently, the critical buckling loads and corresponding buckling mode shapes of the GPLRC shells are obtained. The effects of various factors, including the geometric parameters, boundary conditions and material properties on the torsional buckling behaviors are investigated and discussed in detail.