We present the first Hubble diagram of superluminous supernovae (SLSNe) out to a redshift of two, together with constraints on the matter density, Omega(M), and the dark energy equation-of-state parameter, w(equivalent to 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 data set and improved fitting techniques than previous works. We then solve the SLSN model based on the above standardization via minimization of the chi(2) computed from a covariance matrix that includes statistical and systematic uncertainties. For a spatially flat Lambda cold dark matter (Lambda CDM) cosmological model, we find , with an rms of 0.27 mag for the residuals of the distance moduli. For a w(0)w(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 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 Lambda CDM 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.