A search for the production of three massive vector bosons, $VVZ (V=W, Z)$, in proton-proton collisions at $\sqrt{s} = 13$ TeV is performed using data with an integrated luminosity of $140$ fb$^{-1}$ recorded by the ATLAS detector at the Large Hadron Collider. Events produced in the leptonic final states $WWZ \to \ellν\ellν\ell \ell$ ($\ell=e, μ$), $WZZ \to \ellν\ell\ell \ell\ell$, $ZZZ \to \ell\ell \ell\ell \ell\ell$, and the semileptonic final states $WWZ \to qq \ellν\ell \ell$ and $WZZ \to \ellνqq \ell \ell$, are analysed. The measured cross section for the $pp \rightarrow VVZ$ process is $660^{+93}_{-90}(\text{stat.})^{+88}_{-81}(\text{syst.})$ fb, and the observed (expected) significance is 6.4 (4.7) standard deviations, representing the observation of $VVZ$ production. In addition, the measured cross section for the $pp \rightarrow WWZ$ process is $442 \pm 94 (\text{stat.})^{+60}_{-52}(\text{syst.})$ fb, and the observed (expected) significance is 4.4 (3.6) standard deviations, representing evidence of $WWZ$ production. The measured cross sections are consistent with the Standard Model predictions. Constraints on physics beyond the Standard Model are also derived in the effective field theory framework by setting limits on Wilson coefficients for dimension-8 operators describing anomalous quartic gauge boson couplings.