  We present results from an analysis of all data taken by the BICEP2, Keck
Array and BICEP3 CMB polarization experiments up to and including the 2018
observing season. We add additional Keck Array observations at 220 GHz and
BICEP3 observations at 95 GHz to the previous 95/150/220 GHz data set. The
$Q/U$ maps now reach depths of 2.8, 2.8 and 8.8 $\mu{\mathrm K}_{cmb}$ arcmin
at 95, 150 and 220 GHz respectively over an effective area of $\approx 600$
square degrees at 95 GHz and $\approx 400$ square degrees at 150 & 220 GHz. The
220 GHz maps now achieve a signal-to-noise on polarized dust emission exceeding
that of Planck at 353 GHz. We take auto- and cross-spectra between these maps
and publicly available WMAP and Planck maps at frequencies from 23 to 353 GHz
and evaluate the joint likelihood of the spectra versus a multicomponent model
of lensed-$\Lambda$CDM+$r$+dust+synchrotron+noise. The foreground model has
seven parameters, and no longer requires a prior on the frequency spectral
index of the dust emission taken from measurements on other regions of the sky.
This model is an adequate description of the data at the current noise levels.
The likelihood analysis yields the constraint $r_{0.05}<0.036$ at 95%
confidence. Running maximum likelihood search on simulations we obtain unbiased
results and find that $\sigma(r)=0.009$. These are the strongest constraints to
date on primordial gravitational waves.
