This paper presents an efficient linear system reduction method that computes approximations to the controllability and observability Gramians of large sparse power system descriptor models. The method exploits the fact that a Lyapunov equation solution can be decomposed into low-rank Choleski factors, which are determined by the Alternating Direction Implicit (ADI) method. Advantages of the method are that it operates on the sparse descriptor matrices and does not require the computation of spectral projections onto deflating subspaces of finite eigenvalues, which are needed by other techniques applied to descriptor models. The Gramians, which are never explicitly formed, are used to compute reduced-order statespace models for large-scale systems. Extensions of the method's application to algebraic Riccati equation computation are also considered. Numerical results for small-signal stability power system descriptor models show that the method is efficient for large-scale systems reduced-order model (ROM) computation.