This paper proposes a novel Multivariate Quotient-Difference (MQD) method to obtain the approximate analytical solution for AC power flow equations. Therefore

in the online environment

the power flow solutions covering different operating conditions can be directly obtained by plugging values into multiple symbolic variables

such that the power injections and consumptions of selected buses or areas can be independently adjusted. This method first derives a power flow solution through a Multivariate Power Series (MPS). Next

the MQD method is applied to transform the obtained MPS to a Multivariate Padé Approximants (MPA) to expand the Radius of Convergence (ROC)

so that the accuracy of the derived analytical solution can be significantly increased. In addition

the hypersurface of the voltage stability boundary can be identified by an analytical formula obtained from the coefficients of MPA. This direct method for power flow solutions and voltage stability boundaries is fast for many online applications

since such analytical solutions can be derived offline and evaluated online by only plugging values into the symbolic variables according to the actual operating conditions. The proposed method is validated in detail on New England 39-bus and IEEE 118-bus systems with independent load variations in multi-regions.