Composite materials, such as metal- and polymer-matrix composites, exhibit both elastic (conservative) and viscous (dissipative) features while deforming under com- plex loading histories. Modeling this viscoelastic behavior remains a challenging task, despite the significant advances that have been made in modeling the response of composites with either purely elastic or purely viscous constituents. Even when the phases of the composite are characterized by a single relaxation time (i.e., they ex- hibit a “short term” memory), the viscoelastic interactions between the phases re- sults in a generally continuous spectrum of relaxation times at the macroscopic scale. These “long-term” memory effects manifest themselves as hereditary integrals and contribute significantly to the overall viscoelastic response of the composite in the transient regime (Suquet, 1985). Thus, incorporating these long-term memory effects is crucial for generating accurate estimates for the effective response of viscoelastic composites. In this thesis, a variational procedure is used to generate homogenization estimates for the time-incremental constitutive response of strain-hardening elasto-viscoplastic composites. The generalized standard materials framework (GSM) is exploited to give the local constitutive behavior of the phases a variational representation which is amenable to homogenization. Then, new identities involving the cross-covariances between the local force fields are derived which allow, for the first time, both first and second-order statistical information about the local hardening to be incorporated into the analytical homogenization framework. Next, by taking the limit of the effective time-incremental constitutive relations as the discretization tends to zero, a differential-algebraic system (DAS) of equations characterizing the time-dependent effective constitutive behavior is derived. As a consequence of taking the limit, a long-standing issue of multiple solutions for the LCC properties determining the ef- fective constitutive response is partially resolved. Moreover, the DAS is characterized by primary variables that are physical quantities and are continuous in time, such as the first and second moments of the local fields and their cross-covariances. This novel presentation makes the formulation considerably simpler to implement numeri- cally than the LCC-based formulations, which use fictitious viscosities that can have jump discontinuities, making the time-integration procedure more challenging. It is found that for linear viscoelastic isotropic composites, the DAS recovers the exact results of the correspondence principle for incompressible phases and com- pressible phases under purely hydrostatic loading, while for compressible composites under arbitrary loading the estimates generate fairly accurate estimates (within 1% for monotonic loading and 10% for cyclic loading). The DAS is then used to provide estimates for the effective creep response of nonlinear viscoelastic single-crystals with elastic inclusions. It is found that the presence of elastic (or rigid) inclusions, their morphology and volume fraction have a significant effect on the effective creep rate of the crystal. In particular, the “long-memory” effect manifests in these systems as a transient effective secondary creep-rate (a feature not capture by the Maxwell approximation), and depends strongly on the crystal elasticity and inclusion volume fraction and morphology. Lastly, the model is used to study the interplay between crystallographic anisotropy and morphological anisotropy on the hardening behavior in precipitation strengthened crystals. It is found that precipitates provide a two-pronged strengthening effect on the crystals. On the one hand, they induce a composite effect which increases the macroscopic yield strength, while on the other hand, they promote slip-activity leading to a reduction of the overall anisotropy.