well the question says it all
I have some contour and I want to get the best quadrangular out of it
You have to find an appropriate value for epsilon by using the method of nested intervals, since the number of points increases with decreasing epsion. It may be possible, that the value of 4 corner points cannot be reached though, since the number of points may jump from 3 to 5 at a certain value of epsilon.
If you want to exclude this case, you might need to implement the Ramer–Douglas–Peucker algorithm yourself and modify it, so a given number of points will be returned. (By the way, you're asking for the 'best' 4 corners. You have to specify what best means. approxPolyDp() does not claim an optimal solution!)
Other than that, I see no way to force approxPolyDP() to return 4 points.