In the last century, matter was confirmed to be made up from molecules which consist of two atoms or more. The atom itself consists of a nucleus made of protons and neutrons, and electrons "circling'' around the nucleus. The number of electrons or protons distinguish different elements. Later on, protons and neutrons were found not to be elementary particles but rather composite particles. The question turned then to be what are protons and neutrons made of and this is the focus of elementary particle physics. According to the standard model, protons and neutrons are made up of quarks and gluons. The theory that describes quarks and gluons is called quantum chromodynamics (QCD). According to this theory, quarks and gluons can not be detected freely; they appear only inside hadrons but are never observed freely (confinement). However, at high temperatures and/or densities a transition may happen where quarks and gluons do not exist in bound states (hadrons) anymore but rather exist freely (the asymptotic freedom). This phase of the nuclear matter is known as the quark-gluon plasma (QGP).To learn more about the QCD phase diagram, mainly the confinement and de-confinement transition, many different experiments have been performed from fixed target experiments to high-energy heavy-ion collisions in almost three decades. The discovery of QGP came from ultrarelativistic heavy-ion collision (URHIC) experiments. By ultrarelativistic heavy-ion collisions, we mean heavy ions like gold or lead that have been accelerated to speeds which are close to the speed of light (the ion momentum is much larger than its rest mass). Nowadays, ultrarelativistic heavy-ion collision experiments at the Relativistic Heavy Ion Collider (RHIC) and Large Hadron Collider (LHC) are being used to create and study the quark-gluon plasma. From the early days after confirming the existence of the QGP, relativistic hydrodynamics has been used to describe the hadron spectra and collective flow seen in these experiments and has been quite successful. Since then, different approaches have been developed to model the physics of the QGP. The first approach used was ideal hydrodynamics where the QGP is assumed to behave like a perfect fluid with no viscosity. However, improvements in both the experimental and theoretical sides demonstrated the importance of including dissipative (viscous) effects in QGP modeling. The resulting relativistic viscous hydrodynamics models have been quite successful in describing the data. Despite this success, studies found that the QGP generated in URHICs is a highly momentum-space anisotropic plasma which means that viscous hydrodynamics will break down in some situations. To take this into account, anisotropic hydrodynamics (aHydro) was developed. In aHydro, one includes the momentum-space anisotropies in the distribution function at leading-order, whereas viscous hydrodynamics is expanded around the isotropic distribution function as the leading term and the viscous effects are included as correction terms. In this study, we present a new method for imposing a realistic equation of state in anisotropic hydrodynamics which is called quasiparticle anisotropic hydrodynamics (aHydroQP). In this method, we introduce a single finite-temperature quasiparticle mass which is fit to QCD lattice data. By taking moments of the Boltzmann equation assuming an anisotropic distribution function, we obtain a set of coupled partial differential equations which can be used to describe the 3+1d spacetime evolution of the QGP. Due to the numerical difficulties and the need to understand this new method more, instead of considering the 3+1d case immediately, we begin by studying two simpler cases. First, we specialize to the case of a 0+1d system undergoing boost-invariant Bjorken expansion and compare with the standard method of imposing the equation of state in anisotropic hydrodynamics (aHydro). We find practically no differences between the two methods results for the temperature evolution and the scaled energy density. When we compare the pressure anisotropy, we see only small differences, however, we find significant differences in the evolution of the bulk pressure correction. Second, we present the results in azimuthally-symmetric boost-invariant (1+1d) systems and compare the quasiparticle model with the standard aHydro model and second order viscous hydrodynamics. We compare the three methods' predictions for the primordial particle spectra, total number of charged particles, and average transverse momentum for various values of the shear viscosity to entropy density ratio. We show that they agree well for small shear viscosity to entropy density ratio, but show clear differences at large values of shear viscosity to entropy density ratio. Third, and most importantly, we present the phenomenological predictions of 3+1d quasiparticle anisotropic hydrodynamics compared with LHC 2.76 TeV Pb-Pb collisions. We present comparisons of charged-hadron multiplicity, identified-particle spectra, identified-particle average transverse momentum, charged-particle elliptic flow, identified-particle elliptic flow, elliptic flow as a function of pseudorapidity, and HBT radii. We find good agreement when compared with ALICE data. Looking to the future, we plan to include next-leading-order anisotropic hydrodynamics corrections by including the off-diagonal terms of the anisotropy tensor in quasiparticle anisotropic hydrodynamics. However, since this will be very hard and numerically intense, we consider first next-leading-order anisotropic hydrodynamics using the standard method for imposing the equation of state. To do so, we Taylor-expand assuming small off-diagonal terms to make the formalism easier and numerically tractable. Then, by taking moments of the Boltzmann equation, we find the dynamical equations needed to model the full 3+1d system. In this part of the work, we present only the theory setup and leave the numerical analysis for a future work.