In this paper, we investigate interpolatory model order reduction for large-scale bilinear descriptor systems. Recently, it was shown in [14] for linear descriptor systems that directly extending the standard rational interpolation conditions used in ℋ₂ optimal model reduction to descriptor systems in general yields an unbounded error in the ℋ₂-norm. This is due to the possible mismatch of the polynomial part of the original and reduced-order systems. This conclusion also holds for nonlinear systems as well. In this paper, we deal with bilinear descriptor systems and aim to pay attention to the polynomial part of the bilinear descriptor system along with interpolation. To this end, we have shown in [12] how to determine the polynomial part of each subsystem of the bilinear descriptor system explicitly, by assuming special structures of the system matrices. Considering the same structured bilinear descriptor systems, in this paper we first show how to achieve multipoint interpolation of the underlying Volterra series of bilinear descriptor systems while retaining the polynomial part of each subsystem of the bilinear system. Then, we extend the interpolation based first-order necessary conditions for ℋ₂ optimality to bilinear descriptor systems and propose an iterative scheme to obtain an ℋ₂ optimal reduced-order system. By mean of two numerical examples, we demonstrate the efficiency of the proposed model-order reduction technique and compare it with reduced bilinear systems obtained by using linear IRKA.