Today, more and more substations are created and reconstructed to satisfy the growing electricity demands for both industry and residence. It is always a big concern that the designed substation must guarantee the safety of persons who are in the area of the substation. As a result, the safety metrics (touch voltage, step voltage and grounding resistance), which should be considered at worst case, are supposed to be under the allowable values. To improve the accuracy of calculating safety metrics, at first, it is necessary to have a relatively accurate soil model instead of uniform soil model. Hence, the two-layer soil model is employed in this thesis. The new approximate finite equations with soil parameters (upper-layer resistivity, lower-layer resistivity and upper-layer thickness) are used, which are developed based on traditional infinite expression. The weighted- least-squares regression with new bad data detection method (adaptive weighted function) is applied to fit the measurement data from the Wenner-method. At the end, a developed error analysis method is used to obtain the error (variance) of each parameter. Once the soil parameters are obtained, it is possible to use a developed complex images method to calculate the mutual (self) resistance, which is the induced voltage of a conductor/rod by unit current form another conductor/rod. The basis of the calculation is Green's function between two point current sources, thus, it can be expanded to either the functions between point and line current sources, or the functions between line and line current sources. Finally, the grounding system optimization is implemented with developed three-step optimization strategy using MATLAB solvers. The first step is using "fmincon" solver to optimize the cost function with differentiable constraint equations from IEEE standard. The result of the first step is set as the initial values to the second step, which is using "patternsearch" solver, thus, the non-differentiable and more accurate constraint calculation can be employed. The final step is a backup step using "ga" solver, which is more robust but lager time cost.