Constitutive models are developed for stage I compaction of a two-dimensional array of cylindrical particles under general loading. Densification is assumed to occur by plastic deformation at particle contacts, and yield surfaces are constructed as functions of two state variables using slip-line field and finite element methods. The shape and size of the yield surface are found to depend upon the loading history and the relative density of the array. The effects of non-hydrostatic loading and material hardening upon the pressure-density response of the array of particles are examined by the finite element method. The numerical predictions of the pressure-density response are compared with experimental results for hydrostatic and closed-die compaction of plasticene cylinders.