In this contribution, R-method based approximation of continuous systems to reduced order models (ROMs) is presented by utilizing the grey wolf optimization (GWO) algorithm. The approximation is done by minimizing the errors between time moments (TMs) and Markov parameters (MPs) of the higher-order system (HOS) and desired ROM. So, the TMs and MPs of HOS and ROM are utilized to frame the objective function. The weights associated with objective function are determined using R-method. These weights are further utilized to convert the multi-objective problem into single-objective problem. In objective function, the normalized errors are minimized using GWO algorithm to obtain desired ROM. To ensure the steady-state matching between HOS and its ROM, first TMs of HOS and desired ROM are matched. The stability of obtained ROM is ensured by Hurwitz stability criterion. The superiority of proposed method is determined with the help of two test systems. The results of proposed technique are compared with results of other already obtained ROM available in the literature. For comparative analysis, tabulated values are given for both test cases by considering time domain specifications. These time domain specifications are rise time, settling time, overshoot, undershoot, peak and peak time. The error indices are also provided to validate the proposed method. The step responses, impulse responses, Bode plots and Nyquist plots of the system and ROMs are also presented. The provided results prove the efficacy and effectiveness of the proposed method.