We study the three-dimensional (3D) compact U(1) lattice gauge theory coupled with $N$-flavor Higgs fields by means of the Monte Carlo simulations. This model is relevant to multi-component superconductors, antiferromagnetic spin systems in easy plane, inflational cosmology, etc. It is known that there is no phase transition in the N=1 model. For N=2, we found that the system has a second-order phase transition line $\tilde{c}_1(c_2)$ in the $c_2$(gauge coupling)$-c_1$(Higgs coupling) plane, which separates the confinement phase and the Higgs phase. Numerical results suggest that the phase transition belongs to the universality class of the 3D XY model as the previous works by Babaev et al. and Smiseth et al. suggested. For N=3, we found that there exists a critical line similar to that in the N=2 model, but the critical line is separated into two parts; one for $c_2 < c_{2{\rm tc}}=2.4\pm 0.1$ with first-order transitions, and the other for $ c_{2{\rm tc}} < c_2$ with second-order transitions, indicating the existence of a tricritical point. We verified that similar phase diagram appears for the N=4 and N=5 systems. We also studied the case of anistropic Higgs coupling in the N=3 model and found that there appear two second-order phase transitions or a single second-order transition and a crossover depending on the values of the anisotropic Higgs couplings. This result indicates that an "enhancement" of phase transition occurs when multiple phase transitions coincide at a certain point in the parameter space.