The layer-stacking method can provide three-dimensional conformal dose distributions to the target based on a passive scattering method using mini-spread-out Bragg peak (SOBP). The purpose of this work is to demonstrate the effectiveness of a new weight optimization algorithm that can enhance the robustness of dose distributions against layer depth variation in layer-stacking proton beam therapy. In the robustness algorithm, the upper limit of the layer's weight was adapted to the conventional algorithm and varied for 620 weight set evaluations. The optimal weight set was selected by using an analytical objective function based on Gaussian function with σ = 3 mm for WED variation. Then, we evaluated the stabilities of the one-dimensional depth dose distribution against WED variation generated by Gaussian samples. Three-dimensional dose distributions in the water phantom were also evaluated using the Monte-Carlo dose calculation. The variation of dose as well as dose volume histograms for the spherical target and the organ at risk (OAR) were evaluated. The robustness algorithm reduced the change of the dose distribution due to the WED variation by a factor of almost 3/4 compared to those with the conventional procedure. The rate of 91.8% in total samples was maintained within 5% change of the maximum dose, compared with the rate of 64.9% in the conventional algorithm. In the MC calculation, the high dose-volume in the OAR was reduced around the lateral penumbra and distal falloff region by the robustness algorithm. The stability of depth dose distributions was enhanced under the WED variation, compared to the conventional algorithm. This robust algorithm in layer-stacking proton therapy may be useful for treatment in which the sharpness of the distal falloff along the depth distribution needs to be maintained to spare the organ at risk and keep the dose coverage for the target tumor.