This work focuses on the model order reduction problem for bilinear control systems with nonzero initial conditions. Based on the Volterra series analysis, the system response can be decomposed into three parts. The first two parts are the zero input response and zero initial condition response of the system. The third part describes the response which couples the effect of the nonzero initial condition and the nonzero input. The system corresponding to the third part is a bilinear control system with a special time-varying input coefficient matrix. We show that such a system is equivalent to a time-invariant bilinear control system, and conventional model reduction methods can be applied to reduce it. We propose to reduce each of the component responses independently and then combine them to approximate the full system response. This method is of high flexibility and shows promising results.