This paper is concerned with the shakedown analysis of axisymmetric shells. A new upper-bound nonlinear programming formulation is proposed and solved by an iterative procedure; the solution converges monotonically to an upper bound of the real shakedown load. Moreover, the proposed method is more efficient than the searching procedure in linear programming formulations. A lower-bound linear programming formulation, which is based on Melan's theorem and the stress functions method, is also discussed.