In this work, we have carried out lattice simulations of $(2+1)$-flavor QCD using highly improved staggered quarks at the physical pion mass on $32^3 \times 8$ and $48^3 \times 12$ lattices, with magnetic field strengths ranging up to 0.8 GeV$^2$ and nonzero baryon chemical potentials employing the Taylor expansion framework. We present lattice QCD continuum estimate results, along with the magnetized hadron resonance and ideal gas comparisons, for the leading-order Taylor expansion coefficients for bulk thermodynamic quantities such as pressure, number density, energy density, and entropy density, focusing on the significant impact of strong magnetic fields.