I will present a theory that predicts the viscosity of bimodal and polydisperse colloidal suspensions. Until recently, much work in this area has focused on non-colloidal suspensions where the role of non-hydrodynamic interactions is limited. However, for colloidal dispersions both interparticle forces between pairs of particles and many-particle effects such as depletion forces can have a significant effect on the rheology. As hydrodynamic interactions are also important for colloidal systems, a theoretical description that includes both kinds of interactions is required. An extension of the integral equation theory of Lionberger and Russel (J. Chem Phys. 106 402-416 (1997)) to multicomponent systems accounts for the contribution of thermodynamic interactions to the viscosity of dispersions. Introduction of small particles into a system of larger particles causes depletion forces between the large particles that increase the viscosity, while replacing large particles with an equal volume fraction of small particles increases the free volume available to the system and decreases the viscosity. The integral equations model both of these effects in concentrated suspensions. For a bimodal mixture they predict the dependence of the viscosity on size ratio, composition and total volume fraction. At high volume fraction they locate the minimum of viscosity at composition rich in large particles in reasonable agreement with those measured by Kaler et. al.(Langmuir 8 2382-2389 (1992)). Polydispersity is modeled by a small number of components whose sizes and weights are chosen to match the moments of the size distribution. The predicted decrease in viscosity with polydispersity is in agreement with the Brownian Dynamics simulations of Rastogi et. al. (J. Chem. Phys. 104 9249-9258 (1996)). The hydrodynamic interaction functions that describe the relative motion of a pair of particles are well known for two isolated particles of differing size. We determine the same functions for a pair of particles in a concentrated suspension by distributing the particles according an equilibrium distribution and evaluating the many-body hydrodynamic interactions according to the Stokesian Dynamics simulation method of Brady. For each pair of particles the apparent pair mobility for their separation is computed. This averaging is then repeated over many configurations. At large separations particles behave as if they are in an effective medium with a viscosity dependent on the overall volume fraction of the suspension. From the simulation results, we construct an empirical approximation for these functions based on the three-particle distribution functions. At low volume fractions the addition of hydrodynamic interactions qualitatively changes a viscosity maximum into a viscosity minimum, while at high volume fractions the hydrodynamic interactions simply amplify the trends observed in their absence.