This paper focuses on the systematic optimal tuning of a power system stabilizer (PSS), which can improve the system damping performance immediately following a large disturbance. As the PSS consists of both linear parameters, such as the gain and time constant, and non-smooth nonlinear parameters, such as saturation limits of the PSS, two methods are applied to achieve optimal tuning of all parameters. One is to use the optimization technique based on the Hessian matrix estimated by the feed-forward neural network (FFNN), which identifies the first-order derivatives obtained by the trajectory sensitivities, for the nonlinear parameters. The other is to use an eigenvalue analysis for the linear parameters. The performance of parameters optimized by the proposed method is evaluated by time-domain simulation in both a single-machine infinite bus (SMIB) system and a multi-machine power system (MMPS).