In standard Multidimensional Scaling (MDS) one is concerned with finding a low-dimensional representation of a set of n objects, so that pairwise dissimilarities among the original objects are realized as distances in the embedded space with minimum error. We propose an MDS algorithm that, in addition to minimizing a usual Stress function, can accommodate additional optimization criteria, as well as side constraints associated with the underlying visualization task. We present an application in which we attempt to minimize a secondary objective funcion: the cluster membership discrepancy between a given cluster structure in the original data and the resulting cluster structure in the low-dimensional embedding. Preliminary computational experiments show that the algorithm is able to find MDS embeddings that preserve the original cluster structure while incurring a relatively small increase in Stress, as compared to standard MDS. Finally, we discuss a few properties of the algorithm that make it an interesting choice for Big Data visualization.