Freesteel Blog » Making a molehill out of a circumscribed triangle using numerical minimization algorithms

Making a molehill out of a circumscribed triangle using numerical minimization algorithms

Monday, July 13th, 2015 at 7:05 pm Written by:

Just how bad are these minimization algorithms

I’ve been caught short trusting these classical numerical recipes, packaged up in scipy.optimize.minimize and failing to get an adequate result in the case of the triangular machine tool calibration and its nine unknowns.

So here’s a simple example for finding the circumcircle radius of a triangle whose sides have lengths a, b, c.

The official equation for computing the circumcircle radius is:

def circumcircleradius(a, b, c):
    s = (a+b+c)/2
    areasq = s*(s-a)*(s-b)*(s-c)
    return a*b*c/(4*sqrt(areasq))

To do this using the numerical minimization methods, we consider the triangle having the corner points (0,0), (a,0), (px,py) where (px,py) is the vector we need to solve for to make the other two sides equal in length to values b and c.

def fun1(X):
    px, py = X
    berr = sqrt(px**2 + py**2) - b
    cerr = sqrt((a-px)**2 + py**2) - c
    return berr**2 + cerr**2

Xstartvalue = (a/2, b)
res1 = scipy.optimize.minimize(fun=fun1, x0=Xstartvalue, method='Powell')
px, py = res1.x

That gives us the positions of the three corners of the triangle. Now we need the centre of the circumcircle (cx,cy):

def fun2(X):
    cx, cy = X
    odist = sqrt(cx**2 + cy**2)
    adist = sqrt((cx-a)**2 + cy**2)
    pdist = sqrt((cx-px)**2 + (cy-py)**2)
    # variance formula that gives zero when the three numbers are the same
    var = 3*(odist**2 + adist**2 + pdist**2) - (odist + adist + pdist)**2
    return var

XCstartvalue = (a/2, -b)
res2 = scipy.optimize.minimize(fun=fun2, x0=XCstartvalue, method='Powell')
cx, cy = res2.x

Finally calculate the radius from the distance between the centre and the corner point at (0,0):

rad = sqrt(cx**2 + cy**2)

Now we can test this with a simple right angled triangle where the radius will be half the length of the hypotenuse:

a, b, c = 5, 4, 3
print(numericalcircumcircleradius(a, b, c), circumcircleradius(a, b, c))
>>> 2.5000001856192484 2.5

You can try different numbers, like:

a, b, c = 8, 6, 5
>>> 4.005002372412261 4.005009394574071

Until there’s a failure:

a, b, c = 5, 6, 8
>>> 291816.3453198005 4.005009394574071

Printing the full numerical results:

res1:
  message: 'Optimization terminated successfully.' 
        x: array([-0.30005037, 5.99246565]) 
      nit: 41 
      fun: 8.5873826055365998e-10 

res2:
  message: 'Optimization terminated successfully.' 
        x: array([-232333.60246079, -176572.58157313])
      nit: 121 
      fun: 27.671630859375

The problem is the XCstartvalue starting point is on the wrong side of some boundary where the minimization function thinks it’s a better idea to send it zooming off towards negative infinity from which all three points appear so close together that the distance to each is nearly the same.

Begin with XCstartvalue = (a/2, b) and it’s fine.

If you try to contain the value within a bounding box and use one of the methods which allows one to be set, you can have it completely ignore the limit, or you can sometimes get:

rb = a+b+c
bnds = [(-rb,rb), (-rb,rb)]
res2 = scipy.optimize.minimize(fun=fun2, x0=(a/2, -b), bounds=bnds, method='Nelder-Mead')
>>> message: 'Maximum number of function evaluations has been exceeded.'

Another common attempt to bodge an answer that stays within the limits is to add an upside down table function that tapers upwards outside the region where you know the answer must be:

rb = a+b+c
def fun2(X):
    cx, cy = X
    odist = sqrt(cx**2 + cy**2)
    ...
    if odist > rb:
        var += (odist-rb)**2
    return var

Then:

a, b, c = 5, 6, 8
>>> 19.10269590246688 4.005009394574071

which means it’s found a false minima right up near the edge of the table in the taper zone, since rb = 19.

Using a less smooth function:

def fun2(X):
    ...
    if odist > rb:
        var += abs(odist-rb)
    return var

produces:

a, b, c = 5, 6, 8
>>> 18.99999999998829 4.005009394574071

So that doesn’t work at all. Basically, if within your region there are cells where the solution is going to iterate out towards infinity, then trying to add this table function or force any kind of boundary will merely push such points up against the boundary wrinkle and solve nothing.

What a mess I’ve made from such a trivial example.

These numerical methods solutions are far from the panacea I thought they were going to be.

Leave a comment

XHTML: You can use these tags: <a href="" title=""> <abbr title=""> <acronym title=""> <blockquote cite=""> <code> <em> <strong>