We present the first Hubble diagram of superluminous supernovae (SLSNe) out to a redshift of two, together with constraints on the matter density, ensuremathOmega_M, and the dark energy equation-of-state parameter, w(ensuremathequivp/ensuremathrho). 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 data set and improved fitting techniques than previous works. We then solve the SLSN model based on the above standardization via minimization of the ensuremathchi^2 computed from a covariance matrix that includes statistical and systematic uncertainties. For a spatially flat ensuremathLambda cold dark matter (ensuremathLambdaCDM) cosmological model, we find Omega rm M=0.38^+0.24-0.19, with an rms of 0.27 mag for the residuals of the distance moduli. For a w_0w_aCDM 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 per cent. 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 ensuremathLambdaCDM model with the same precision (considering only statistical uncertainties) as current surveys that use Type Ia supernovae, while providing a factor of 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.