Based on Melan’s theorem, an improved numerical solution procedure for evaluating shakedown loads by non-linear superposition method is presented, and the relationship between the classical non-linear programming formulation of shakedown problem and the numerical method is disclosed. The stress term in classical optimization problem is replaced by the element bearing-ratio (EBR) in the procedure, and series of residual EBR fields can be generated by the D-value of the elastic-plastic EBR fields and the elastic EBR fields at every incremental loading step. The shakedown load is determined by performing the incremental non-linear static analysis when the yield criterion is arrived either by the elastic-plastic EBR fields or residual EBR fields. By introducing the EBR, the proposed procedure can be easily used to those complex structures with multi-material and complicated configuration. The procedure is described in detail and some numerical results, that show the efficiency and effectiveness of the proposed method, are reported and discussed.