We present a simulation-based cosmological analysis using a combination of Gaussian and non-Gaussian statistics of the weak lensing mass (convergence) maps from the first three years of the Dark Energy Survey. We implement the following: (1) second and third moments; (2) wavelet phase harmonics; (3) the scattering transform. Our analysis is fully based on simulations, spans a space of seven w Cold Dark Matter (wCDM) cosmological parameters, and forward models the most relevant sources of systematics inherent in the data: masks, noise variations, clustering of the sources, intrinsic alignments, and shear and redshift calibration. We implement a neural network compression of the summary statistics, and we estimate the parameter posteriors using a simulation-based inference approach. Including and combining different non-Gaussian statistics is a powerful tool that strongly improves constraints over Gaussian statistics (in our case, the second moments); in particular, the figure of merit (S8,Ωm) is improved by 70% (ΛCDM) and 90% (wCDM). When all the summary statistics are combined, we achieve a 2% constraint on the amplitude of fluctuations parameter S8≡σ8(Ωm/0.3)0.5, obtaining S8=0.794±0.017 (ΛCDM) and S8=0.817±0.021 (wCDM), and a ∼10% constraint on Ωm, obtaining Ωm=0.259±0.025 (ΛCDM) and Ωm=0.273±0.029 (wCDM). In the context of the wCDM scenario, these statistics also strengthen the constraints on the parameter w, obtaining w<-0.72. The constraints from different statistics are shown to be internally consistent (with a p-value>0.1 for all combinations of statistics examined). We compare our results to other weak lensing results from the first three years of the Dark Energy Survey data, finding good consistency; we also compare with results from external datasets, such as Planck constraints from the cosmic microwave background, finding statistical agreement, with discrepancies no greater than <2.2σ.