We present a numerical method for the computation of shakedown loads of structures under thermo-mechanical loading accounting for limited kinematical hardening. The method is based on the lower bound approach by Melan extended to the hardening using a two-surface model. Both the yield and the bounding surface are defined by the von Mises condition. Melan’s shakedown theorem leads to a nonlinear convex optimization problem. This is solved by the interior-point algorithm ipsa recently developed by the authors. In this paper, theoretical and numerical aspects will be presented as well as numerical examples from mechanical engineering.