We study strongly correlated superconductivity in magic-angle twisted bilayer graphene (MATBG) using a variational Gutzwiller wavefunction |Ψ_G⟩ = ∏_R P̂_R |Φ₀⟩, where the Gutzwiller projector P̂_R is allowed to break charge U(1) symmetry to accommodate superconducting (SC) order. The ground state energy is evaluated via the Gutzwiller Approximation applied to an 8-band model consisting of correlated f-orbitals and uncorrelated c-orbitals, with interactions including onsite Coulomb repulsion U, phonon-mediated anti-Hund's coupling Ĥ_(J_A), and intra-orbital Hund's coupling Ĥ_(J_H). At filling ν = 2.5, we map out the phase diagram as a function of U and J_A, and reveal a strongly correlated SC (SC-SC) phase dominates at large U, wherethe strong on-site interaction U strongly suppress the f-orbital charge fluctuations while maintaining finite pairing order and a sizeable quasiparticle weight Z, distinguishing it from a conventional Mott insulator. For a range of J_(rm A), SC-SC transitions to FL as U decreases, until the weakly correlated BCS-like SC (BCS-SC) re-enters as U → 0. We further identify a novel small Fermi liquid (sFL) state with effective Fermi surface formed by c-orbitals, which is essentially different with the normal Fermi liquid. Interestingly, in the intermediate- (U ≲ 40 meV) and large-U (U ≳ 40 meV) regimes, the conventional FL and the sFL are the lowest-energy normal phases, respectively, potentially serve as the parent states of the SC-SC phase. These results illuminate the interplay between strong correlations and unconventional pairing in MATBG, and establish a versatile Gutzwiller framework applicable to other strongly correlated superconductors.