This　paper　presents　an　efficient　operational　matrix　method　for　two　dimensional　nonlinear　variable　order　fractional　optimal　control　problems　(2D-NVOFOCP).　These　problems　include　the　nonlinear　variable　order　fractional　dynamical　systems　(NVOFDS)　described　by　partial　differential　equations　such　as　the　diffusion-wave,　convection-diffusion-wave　and　Klein-Gordon　equations.　The　variable　order　fractional　derivative　is　defined　in　the　Caputo　type.　The　proposed　hybrid　method　is　based　on　the　transcendental　Bernstein　series　(TBS)　and　the　generalized　shifted　Chebyshev　polynomials　(GSCP).　The　new　operational　matrices　of　derivatives　are　generated　for　the　mentioned　polynomials.　The　state　and　control　functions　are　expressed　by　the　TBS　and　GSCP　with　free　coefficients　and　control　parameters.　These　expansions　are　substituted　in　the　performance　index　and　the　resulting　operational　matrices　are　employed　to　extract　algebraic　equations　from　the　approximated　NVOFDS.　The　constrained　extremum　is　obtained　by　coupling　the　algebraic　constraints　from　the　dynamical　system　and　the　initial　and　boundary　conditions　with　the　algebraic　equation　extracted　from　the　performance　index　by　means　of　a　set　of　unknown　Lagrange　multipliers.　The　convergence　analysis　is　discussed　and　several　numerical　experiments　illustrate　the　efficiency　and　accuracy　of　the　proposed　method.　©　2020　Chinese　Automatic　Control　Society　and　John　Wiley　&　Sons　Australia,　Ltd