Homogeneous-in-the-state bilinear systems, appended by an additive disturbance, appear both from the discretization of some partial differential equations and from the bilinearization of certain nonlinear systems. They often have large state vectors that can be cumbersome for simulation and control system design. Our aim is to define a method, invariant to time transformations, for finding a reduced-order model with similar disturbance-output characteristics to those of the plant for all admissible input sequences. The inputs considered satisfy simple upper and lower bound constraints, representing saturating actuators. The approximation is based on a model truncation approach and a condition for the existence of such an approximation is given in terms of the feasibility of a set of linear matrix inequalities. A novelty of our work is in the definition of a new Gramian for this class of systems. Explicit error bounds on the scheme are included. The paper concludes with a demonstration of the proposed approach to the model reduction of a solar collector plant. (C) 2011 Elsevier Ltd. All rights reserved.