Approximating polynomial under monotonicity constraint

Question

Given a bivariate data set, how does one obtain the most accurate least-squares polynomial approximation under the constraint that the resulting polynomial is non-decreasing. Similarly, how does one proceed with a trivariate data set with the same constraint applying to the resulting bivariate polynomial.