Microbubbles coated by visco-elastic shells are important for ultrasound diagnosis using contrast agents, and dynamics of single coated bubble has been investigated in literature. However, although a high number of contrast agents are used in practical situations, there has long been an absense of nonlinear acoustic theory for multiple coated bubbles, except for our recent work (Kikuchi & Kanagawa, Jpn. J. Appl. Phys., 60 (2021), SDDD14), under several assumptions to be excluded. Aiming for generalization, in this study, we theoretically investigate weakly nonlinear propagation of ultrasound in liquid containing multiple bubbles coated by a visco-elastic shell with compressibility. Leveraging the method of multiple scales, both the Korteweg-de Vries-Burgers (KdVB) equation for a low-frequency long wave and nonlinear Schrödinger (NLS) equation for a high-frequency short wave, are derived from the volumetric averaged equations for bubbly liquids based on a two-fluid model and the up-to-date model for single coated bubble with shell compressibility. Neglected factors in our previous paper, i.e., compressibility of the shell and liquid, drag force acting on bubbles, bubble translation, and thermal conduction, are incorporated in the present KdVB and NLS equations; the proposed model will be regarded as a generic physico-mathematical model. Results show that shell compressibility attenuated ultrasound strongly and decreased nonlinearity of ultrasound. Finally, we compared the magnitudes of six dissipation factors (shell compressibility, shell viscosity, liquid compressibility, liquid viscosity, thermal effect, and drag force) for five typical ultrasound contrast agents, and a similar tendency between KdVB and NLS equations was revealed.