  In this paper, we consider model order reduction for bilinear systems with
non-zero initial conditions. We discuss choices of Gramians for both the
homogeneous and the inhomogeneous parts of the system individually and prove
how these Gramians characterize the respective dominant subspaces of each of
the two subsystems. Proposing different, not necessarily structure preserving,
reduced order methods for each subsystem, we establish several strategies to
reduce the dimension of the full system. For all these approaches, error bounds
are shown depending on the truncated Hankel singular values of the subsystems.
Besides the error analysis, stability is discussed. In particular, a focus is
on a new criterion for the homogeneous subsystem guaranteeing the existence of
the associated Gramians and an asymptotically stable realization of the system.
