This paper proposes a numerical solution method for upper bound shakedown analysis of perfectly elasto-plastic thin plates by employing the C1 natural element method.Based on the Koiter's theorem and von Mises yield criterion,the nonlinear mathematical programming formulation for upper bound shakedown analysis of thin plates is established.In this formulation,the trail function of residual displacement increment is approximated by using the C1 shape functions,the plastic incompressibility condition is satisfied by introducing a constant matrix in the objective function,and the time integration is resolved by using the K(o)nig's technique.Meanwhile,the objective function is linearized by distinguishing the non-plastic integral points from the plastic integral points and revising the objective function and associated equality constraints at each iteration.Finally,the upper bound shakedown load multipliers of thin plates are obtained by direct iterative and monotone convergence processes.Several benchmark examples verify the good precision and fast convergence of this proposed method.