Abstract Non-linear stability plays an essential role in the design of reticulated shell structures. In order to improve the buckling strength of aluminum alloy spherical shells with gusset joints, this paper proposes a shape optimization method to maximize the non-linear buckling load using a genetic algorithm. The control points of a cubic spline are chosen as the design variables defining the height of the surface. In the non-linear buckling analysis, a linear combination of the symmetric and antisymmetric buckling modes is used as the imperfection mode, and an initial eccentricity is introduced for each member to consider the disadvantageous effect caused by the interaction between the flexural-torsional buckling and global buckling. It is verified in the numerical examples that the proposed method can increase the non-linear buckling load maintaining the smoothness of the surface, while the conventional approach of minimizing the total strain energy is not effective for improving the buckling strength. Finally, the design tables of the parameters of the optimized cubic spline curve are obtained through the parametric analysis considering various load factors and height-to-span ratios, and are verified to be practical and useful under different height of the section, number of rings and support condition.