An energy functional given by W = integral integral integral d3x[B2/(2mu0) + p(parallel-to)/(GAMMA - 1)] is proposed as a variational principle to determine three-dimensional (3D) magnetohydrodynamic (MHD) equilibria with anisotropic plasma pressure. It is demonstrated that the minimisation of W using an inverse coordinate spectral method reproduces the force balance relations that govern the MHD equilibrium properties of 3D plasmas with p(parallel-to) not-equal p(parpendicular-to) that have nested magnetic flux surfaces. Numerical procedures already developed for the scalar pressure model can be easily extended to the anisotropic pressure model. Specifically, a steepest descent procedure coupled with the application of a preconditioning algorithm to improve the convergence behaviour has been employed to minimise the energy of the system. The numerical generation of 3D torsatron equilibria with highly localised anisotropic pressure distributions attests to the robustness of the method of solution considered.