In this paper, we present a new technique by combining the Taylor series expansion with the Arnoldi method to automatically develop reduced-order models for coupled energy domain nonlinear microelectromechanical devices. An electrostatically actuated fixed-fixed beam structure with squeeze-film damping effect is examined to illustrate the model-order reduction method. Simulation results show that the reduced-order nonlinear models can accurately capture the device dynamic behavior over a much larger range of device deformation than the conventional linearized model. Compared with the fully meshed finite-difference method, the model reduction method provides accurate models using orders of magnitude less computation. The reduced MEMS device models are represented by a small number of differential and algebraic equations and thus can be conveniently inserted into a circuit simulator for fast and efficient system-level simulation.