Continuation of periodic solutions in dynamical systems can be efficiently carried out using the harmonic balance method (HBM). Various tools exist for the stability analysis of HBM solutions of smooth systems, such as the Koopman-Hill method, but many real-world phenomena are modeled more accurately by nonsmooth dynamics. We generalize the Koopman-Hill stability method to constrained multibody systems undergoing planar Coulomb friction. The harmonic balance method has recently been applied successfully to systems subject to nonsmooth Coulomb friction. The system dynamics together with the set-valued friction law, expressed as an algebraic constraint equation, lead to a differential-algebraic equation (DAE) of changing index. We show that the Jacobian of the HBM residual, also known as the Hill matrix, defines a lifted linear time-invariant DAE, in analogy to the ODE case. This DAE can be solved in closed form using a generalized matrix exponential involving Drazin inverses. The resulting monodromy matrix is singular, indicating that only admissible initial conditions can be evolved over a period. Applied to a frictional two-mass oscillator, the generalized Koopman-Hill method reliably identifies Floquet multipliers that align with a reference obtained through event-driven integration. The spectrum of the Hill matrix is fragmented due to the nonsmooth character of the solutions, which makes it difficult to identify Floquet exponents using sorting methods. However, the proposed generalized Koopman-Hill method reliably identifies the correct Floquet exponents.