Aiming to address the limitation of current non-probabilistic reliability methods for slopes, which can only construct an uncertain domain with regular boundaries and a large envelope when using the convex set model to describe parameter uncertainty, this study proposes a novel non-probabilistic reliability analysis method for slopes based on the multi-convex set model. The approximate performance function of the slope is formulated by the quadratic response surface method in conjunction with Latin hypercube sampling. Furthermore, both the traditional interval model and the PCA (principal component analysis)-based interval model are established. By integrating these two models, a multi-convex set model is constructed. The HL-RF (Hasofer-Lind and Rackwitz-Fiessler) iterative algorithm is employed to identify the most probable failure point of the limit state function, while the simplex optimization algorithm is utilized to locate the extreme point. Based on the definition of the non-probabilistic reliability index as a distance ratio, the non-probabilistic reliability of the slope is calculated, and its stability status is assessed accordingly. The feasibility of the proposed method is validated through case studies. Compared with non-probabilistic reliability methods for slopes based on the interval and ellipsoid model, the results obtained by the presented method exhibit greater consistency with those derived from the Monte Carlo method. As the variability and correlation of shear strength parameters increase, the non-probabilistic reliability index of the slope decreases. When applied to slope stability analysis, the judgment outcomes align well with those obtained via various reliability approaches.