Nonlinear balanced truncation is a model order reduction technique that reduces the dimension of nonlinear systems in a manner that accounts for either open- or closed-loop observability and controllability aspects of the system. A computational challenges that has so far prevented its deployment on large-scale systems is that the energy functions required for characterization of controllability and observability are solutions of various high-dimensional Hamilton–Jacobi–(Bellman) equations, which are computationally intractable in high dimensions. This work proposes a unifying and scalable approach to this challenge by considering a Taylor-series-based approximation to solve a class of parametrized Hamilton–Jacobi–Bellman equations that are at the core of nonlinear balancing. The value of a formulation parameter provides either open-loop balancing or a variety of closed-loop balancing options. To solve for the coefficients of Taylor-series approximations to the energy functions, the presented method derives a linear tensor system and heavily utilizes it to numerically solve structured linear systems with billions of unknowns. The strength and scalability of the algorithm is demonstrated on two semi-discretized partial differential equations, namely the Burgers and the Kuramoto–Sivashinsky equations.