  We present the first Hubble diagram of superluminous supernovae (SLSNe) out
to a redshift of two, together with constraints on the matter density,
$\Omega_{\rm M}$, and the dark energy equation-of-state parameter, $w(\equiv
p/\rho)$. We build a sample of 20 cosmologically useful SLSNe~I based on light
curve and spectroscopy quality cuts. We confirm the robustness of the peak
decline SLSN~I standardization relation with a larger dataset and improved
fitting techniques than previous works. We then solve the SLSN model based on
the above standardisation via minimisation of the $\chi^2$ computed from a
covariance matrix which includes statistical and systematic uncertainties. For
a spatially flat $\Lambda$CDM cosmological model, we find $\Omega_{\rm
M}=0.38^{+0.24}_{-0.19}$, with a rms of 0.27 mag for the residuals of the
distance moduli. For a $w_0w_a$CDM cosmological model, the addition of SLSNe~I
to a `baseline' measurement consisting of Planck temperature together with type
Ia supernovae, results in a small improvement in the constraints of $w_0$ and
$w_a$ of 4\%. We present simulations of future surveys with 868 and 492 SLSNe I
(depending on the configuration used) and show that such a sample can deliver
cosmological constraints in a flat $\Lambda$CDM model with the same precision
(considering only statistical uncertainties) as current surveys that use type
Ia supernovae, while providing a factor 2-3 improvement in the precision of the
constraints on the time variation of dark energy, $w_0$ and $w_a$. This paper
represents the proof-of-concept for superluminous supernova cosmology, and
demonstrates they can provide an independent test of cosmology in the
high-redshift ($z>1$) universe.
