There may be one to an infinite number of solutions. A simple example of infinitely many solutions is a case like orthogonal rectangles with the edge of one rectangle being closest and equidistant to the other. If the surfaces intersect in a curve there are also infinitely solutions. There are also many situations where there are multiple discrete solutions and these occur with disturbing frequency in the real world of modeling.
You need to carefully consider what you want to do in these cases.
Typically you will use an iterative techniques find a local solution, and then look for a better one until you are certain no better one exists.
Mathematically, you are searching for an extreme value of a function of four parameters over a bounded domain. The function is differentiable, but the derivatives can be parallel or zero (when control points are coincident.) To handle cases when the closest point or points are on a boundary, you will need code that finds extreme points of functions with two, three and four parameters.
This is a tough numerical problem. If you have not written this kind of code before, read the optimization sections of something like Numerical Recipes in C.
I typically use a three phase approach to this style of problem
First: A fast and simple algorithm that gets a decent answer most of the time (bounded Newton style stuff) together with a reliable way to determine that the fast and simple algorithm failed. Detecting failure is the tricky part. The algorithm can converge and still fail miserably when derivatives are parallel-ish or zero-ish (the difference between precision and accuracy). When these simple algorithms work, and they do work most of time on “reasonable” surfaces, they tend to be fast and accurate.
Second: A slower but better algorithm to try again (conjugate gradient, etc.) and a reliable way to determine if it failed. Again, detecting failure is the tricky part
Third: A guided sampling style algorithm as a last ditch fallback (these are typically much slower than bounded newton stuff because they require lots of iterations.)
I often use course tessellations to get seed values for iterative algorithms.
We need to expose the closest distance functions we have in the opennurbs SDK Rhino plug-in devs have access to, so our customers don’t need to write stuff like this. I’ve added http://mcneel.myjetbrains.com/youtrack/issue/RH-28493 to remind me to get this done.