  Interpolatory projection methods for model reduction of nonparametric linear
dynamical systems have been successfully extended to nonparametric bilinear
dynamical systems. However, this is not the case for parametric bilinear
systems. In this work, we aim to close this gap by providing a natural
extension of interpolatory projections to model reduction of parametric
bilinear dynamical systems. We introduce necessary conditions that the
projection subspaces must satisfy to obtain parametric tangential interpolation
of each subsystem transfer function. These conditions also guarantee that the
parameter sensitivities (Jacobian) of each subsystem transfer function is
matched tangentially by those of the corresponding reduced order model transfer
function. Similarly, we obtain conditions for interpolating the parameter
Hessian of the transfer function by including extra vectors in the projection
subspaces. As in the parametric linear case, the basis construction for
two-sided projections does not require computing the Jacobian or the Hessian.
