By coupling the Boltzmann transport equations of both quarkonium and open heavy quarks, we investigate their dynamical evolution inside the quark-gluon plasma and study quarkonium production in heavy ion collisions. The Boltzmann transport equation of quarkonium is derived from the open quantum system formalism and effective field theory of QCD by assuming quarkonium interacts weakly with the plasma. The dissociation and recombination terms in the Boltzmann equation are calculated in potential nonrelativistic QCD. It is shown that the combination of quarkonium dissociation, recombination, open heavy quark diffusion and energy loss can drive the system of quarkonium and open heavy quarks to detailed balance and kinetic thermalization. By solving the transport equations with initial momenta of quarkonia and heavy quarks sampled from PYTHIA and a hydrodynamic medium, we can calculate the nuclear modification factors of bottomonium and describe the data at both RHIC and LHC energies. The azimuthal angular anisotropy coefficient $v_2$ of $\Upsilon$(1S) in 5.02 TeV peripheral Pb-Pb collisions is also predicted.