A challenge in many applications of non-parametric curve estimation is that the function must satisfy some (lower and/or upper) variable order constraints (for example, a density is constrained to lie between two functions). At the same time the spatially inhomogeneous smoothness of the function is modelled by Besov and Triebel-type smoothness constraints. Donoho and Johnstone (1998) and Delyon and Juditsky (1996) studied minimax rates of convergence for wavelet estimators with thresholding, while Lepski et al. (1997) proposed a variable bandwidth selection for kernel estimators that achieved optimal rates over the scale of Besov spaces. Here we show how to construct estimators that satisfy the variable order constraints and also achieve minimax rates over the appropriate smoothness class. This generalizes results of Dechevsky and MacGibbon (1999) for the case of constant constraints. The parameters of the new constrained estimator (when the constraints are functions) are shown here to depend on the regularity of the constraint functions, except when the lower constraint function is convex and/or the upper constraint function is concave. A preliminary announcement of some of the results of the present work (without proofs) was made in Dechevsky (2007) as a part of a survey on the state of the art and ongoing research in shape-preserving wavelet approximation.