The paper develops a method for model reduction of bilinear control systems. It leans upon the observation that the input-output map of a bilinear system has a particularly simple Fliess series expansion. Subsequently, a model reduction algorithm is formulated such that the coefficients of Fliess series expansion for the original and reduced systems match up to certain predefined sets -nice selections. Algorithms for computing matrix representations of unobservability and reachability spaces complying with a nice selection are provided. Subsequently, they are used for calculating a partial realization of a given input-output map.Prior work To the best of our knowledge, the contribution of the paper is new. Model reduction of bilinear systems is an established topic, without claiming completeness, we mention Bai and Skoogh (2006); Zhang and Lam (2002); Wang and Jiang (2012); Xu et al. (2015); Breiten and Damm (2010); Benner and Breiten (2015); Flagg and Gugercin (2015); Feng and Benner (2007); Lin et al. (2007); Flagg (2012). In particular, moment matching methods for bilinear systems were proposed in Bai and Skoogh (2006); Feng and Benner (2007); Lin et al. (2007); Flagg (2012); Breiten and Damm (2010); Benner and Breiten (2015); Flagg and Gugercin (2015); Wang and Jiang (2012) and for general non-linear systems in Astolfi (2010). The current paper represents another version of moment matching for bilinear systems. In Flagg (2012); Breiten and Damm (2010); Benner and Breiten (2015) the concept of moment matching at some frequencies σ 1 , . . . , σ k was defined. The papers Bai and Skoogh (2006); Feng and Benner (2007); Lin et al. (2007) correspond to moment matching at zero frequency (σ = 0). The cur-arXiv:1605.04414v1 [math.OC]