This article presents an elastoplastic continuum topology optimization method for shakedown under the SIMP‐based framework for the first time, aiming to maximize shakedown load‐carrying capacity under variable repeated loading. In contrast to most elastoplastic topology optimizations which have path‐dependent properties, the shakedown topology optimization only needs loading vertices of the loading domain rather than any complete loading history. Based on Melan's lower bound theorem, the gradient‐based topology optimization framework synthesizing the shakedown analysis and sensitivity analysis is developed to maximize the shakedown multiplier. Path‐independent sensitivity related to the specified density interpolation model is derived analytically for the first time in conjunction with the adjoint method. Besides, a primal‐dual algorithm of shakedown analysis is implemented for computing the shakedown load limit of structures and corresponding adjoint variables needed in the sensitivity analysis. In addition to the density interpolation scheme of material Young's module, the yield strength is also interpolated to overcome numerical difficulties. Several numerical examples demonstrate the effectiveness of the proposed method, which indicates the shakedown load‐carrying performance enhancements.