The conventional portfolio design approach assumes Gaussian return distributions, but this is not accurate in practice. Asymmetric and heavy-tailed return distributions necessitate consideration of higher-order moments such as skewness and kurtosis, in addition to mean and variance. This study proposes a multi-objective approach using a mean-variance-skewness-kurtosis model to construct a diversified portfolio. A parametrized polynomial goal programming (PGP) method is used to optimize the portfolio by maximizing returns and skewness while minimizing variance and kurtosis. Empirical data from the S &P ESG index family is used, and PGP generates multiple portfolios reflecting investors’ preferences for the four moments. To compare between the obtained portfolios, we represent the empirical cumulative distribution of the portfolio returns for all studied values of weights and show how this can be used to assist the inverstor in selecting the best set of weights.