Fermat-Torricelli problem redux

The Fermat-Torricelli problem is one of the most interesting in optimization theory. The classical version of this problem is about finding the point inside a triangle that is the minimum distance to its three vertices. The problem was proposed by Pierre de Fermat and solved by Evangelista Torricelli, who provided actually several solutions to the problem. A more generalized version of the problem tries to find the point \((x^\ast,y^\ast)\) in a triangle with vertices \((x_1,y_1)\), \((x_2,y_2)\), and \((x_3,y_3)\) that minimizes the function \[ F(x,y)=\sum_{i=1}^3 w_i \sqrt{(x_i-x)^2+(y_i-y)^2} \] where \(w_i\) are weights. When these weights are all equal to one, minimizing function \(F(x,y)\) solves the classic Fermat-Torricelli problem. When the weights are different, minimizing \(F(x,y)\) solves a transportation problem that is also known as the Steiner problem, or the railway junction problem. In this case, point \((x^\ast,y^\ast)\) minimizes transportation costs when \(w_i\) are the volume of traffic originating from vertex \(i\).

Several approaches are known that find \((x^\ast,y^\ast)\) numerically and iteratively. But recently a Russian mathematician, Alexei Yur'evich Uteshev from Saint-Petersburg State University, published a paper Analytical Solution to the Generalized Fermat-Torricelli Problem in The American Mathematical Monthly 121(4) April 2014 that provides a nifty analytical solution to the problem. An earlier version of his paper is also available on arXiv.

P1 P2 P3 P* The illustration on the right shows an example. A triangle with three vertices (P1, P2, P3) can be thought of as three cities. We now want to build three highways originating in each city and meet in P*. The thickness of the green lines reflects the capacity each highway segement is expected to handle. As is visible in the illustration, the larger traffic volume originating in P3 "pulls" the point P* where the highways meet closer to P3. Because there is also more traffic originating from P2 than P1, the junction P* is closer to the line segement P2–P3 than to P1–P3.

I got interested in the Fermat-Torricelli problem because it is solves a class of transport optimization problems, also in electric grid design. If the weights in the above problem reflect traffic volume (or electric loads), then \(w_1= v_{12}+v_{13}\), \(w_2= v_{12}+v_{23}\), and \(w_3= v_{13}+v_{23}\). In simple terms, each of the three segments that connect to point P* must handle the sum of the traffic that forks off to the other two points. It is interesting to note that this variation of the problem can have solutions where P* is almost (but not exactly) on one of the line segments that connects the triangle vertices, with the "spur line" branching off perpendicularly. For example, if \(v_{23}\) is very large compared to \(v_{12}+v_{13}\), then \(w_1\) is small and both \(w_2\) and \(w_3\) are relatively large.

Another interesting variation of the Fermat-Torricelli problem is to think about the redundancy that is built into a network that connected all three vertices (the blue lines in the diagram above). If a network is connecting through a central hub, failure of any of the hub connectors (the green lines in the diagram) stops all traffic. While a network consisting of only blues lines is more expensive, it tends to be more robust. For example, if P2–P3 fails, traffic from P2 to P3 can be rerouted through P1. A network that is prone to disruptions would probably favour this then a less costly but "buggy" network. The cost of outage would have to be added to the optimization problem to get the "right" answer. In other words, uncertainty matters when it comes to economical solutions.

Below is a code fragment in the SAS matrix language that reads in a database of triangle coordinates and weights and writes out a database with the Fermat-Torricelli points. The code also handles the cases where an interior solution does not exist. When the weights for one vertex are too large, that vertex may become the optimal point. Intuitively, if most traffic is to and from point P3, it is not necessary to build a highway between P1 and P2.

The code below tries to be efficient about some of the calculations. The distances between points appear as their squares (e.g., r12s) rather than as actual distances. This speeds up computation. The output from the IML code provides the coordinates of the optimal Fermat-Torricelli point \((x^\ast,y^\ast)\) as well as the value of \(F(x^\ast,y^\ast)\). Below the SAS macro are a few sample calculations for testing. Example9 corresponds to the points in the SVG diagram above.

*------------------------------------------------------------+ | Calculation of Generalized Fermat-Torricelli Points | | Following "Analytical Solution for the Generalized | | Fermat-Torricelli Problem" by Alexei Yu. Uteshev | +------------------------------------------------------------+ | (C) 2014 Werner Antweiler, University of British Columbia | +------------------------------------------------------------; %MACRO FERMATPOINT(DBASE,ID,X1,Y1,M1,X2,Y2,M2,X3,Y3,M3); proc iml; use &DBASE; read all var { &id } into id; read all var { &X1 &X2 &X3 } into x; read all var { &Y1 &Y2 &Y3 } into y; read all var { &M1 &M2 &M3 } into m; close &DBASE; n=nrow(X); start euclidean2(x1,y1,x2,y2); return(((x1-x2)**2+(y1-y2)**2)); finish; start cos_angle(r0s,r1s,r2s); return((r1s+r2s-r0s)/(2*sqrt(r1s*r2s))); finish; items={"X_Fermat","Y_Fermat","Fermat_min"}; XY=J(n,3,.); do i=1 to n; r12s=euclidean2(x[i,1],y[i,1],x[i,2],y[i,2]); r23s=euclidean2(x[i,2],y[i,2],x[i,3],y[i,3]); r13s=euclidean2(x[i,3],y[i,3],x[i,1],y[i,1]); cos12=cos_angle(r12s,r23s,r13s); cos23=cos_angle(r23s,r12s,r13s); cos13=cos_angle(r13s,r12s,r23s); u1=m[i,2]**2+m[i,3]**2+2*m[i,2]*m[i,3]*cos23-m[i,1]**2; u2=m[i,1]**2+m[i,3]**2+2*m[i,1]*m[i,3]*cos13-m[i,2]**2; u3=m[i,1]**2+m[i,2]**2+2*m[i,1]*m[i,2]*cos12-m[i,3]**2; if (u1<=0) then do; xy[i,1]=x[i,1]; xy[i,2]=y[i,1]; xy[i,3]=m[i,2]*sqrt(r12s)+m[i,3]*sqrt(r13s); end; else if (u2<=0) then do; xy[i,1]=x[i,2]; xy[i,2]=y[i,2]; xy[i,3]=m[i,1]*sqrt(r12s)+m[i,3]*sqrt(r23s); end; else if (u3<=0) then do; xy[i,1]=x[i,3]; xy[i,2]=y[i,3]; xy[i,3]=m[i,1]*sqrt(r13s)+m[i,2]*sqrt(r23s); end; else do; S=abs(x[i,1]*y[i,2]+x[i,2]*y[i,3]+x[i,3]*y[i,1] -x[i,1]*y[i,3]-x[i,3]*y[i,2]-x[i,2]*y[i,1]); MMM=(m[i,1]+m[i,2]+m[i,3])/2; MMS=(MMM*(MMM-m[i,1])*(MMM-m[i,2])*(MMM-m[i,3])); sigma=2*sqrt(MMS); K1=(r12s+r13s-r23s)*sigma+(m[i,2]**2+m[i,3]**2-m[i,1]**2)*S; K2=(r23s+r12s-r13s)*sigma+(m[i,1]**2+m[i,3]**2-m[i,2]**2)*S; K3=(r13s+r23s-r12s)*sigma+(m[i,1]**2+m[i,2]**2-m[i,3]**2)*S; d=(m[i,1]**2*K1+m[i,2]**2*K2+m[i,3]**2*K3)/(2*sigma); KS=K1*K2*K3/(4*S*sigma*d); xy[i,1]=KS*(x[i,1]/K1+x[i,2]/K2+x[i,3]/K3); xy[i,2]=KS*(y[i,1]/K1+y[i,2]/K2+y[i,3]/K3); xy[i,3]=sqrt(d); end; end; create fermatpoint from XY[r=id c=items]; append from XY[r=id]; close fermatpoint; quit; run; %MEND; options nocenter ls=80; data sample; input id $ x1 y1 x2 y2 x3 y3 m1 m2 m3; datalines; Example1 2 6 1 1 5 1 2 3 4 Example2 2 6 1 1 5 1 3 5 4 Example3 0 0 2 0 -1.41421356 1.41421356 1.5 2 2 Example4 39 57 22 42 42 75 18 41 52 Example5 -0.5 -0.288675 0.5 -0.288675 0 0.57735 1 1 1 Example6 -0.5 -0.288675 0.5 -0.288675 0 0.57735 1 1 1.732 Example7 -0.5 -0.288675 0.5 -0.288675 0 0.57735 1 1 2 Example8 -0.5 -0.288675 0.5 -0.288675 0 0.57735 1 1 0.01 Example9 25 90 120 210 190 30 3 4 5 ; %FERMATPOINT(sample,id,x1,y1,m1,x2,y2,m2,x3,y3,m3); proc print data=fermatpoint; id id; var X_fermat Y_Fermat Fermat_min; run;

With this somewhat esoteric blog written on Christmas eve, I wish everyone a Merry Christmas 2014. May the world find peace and love.

Posted on Wednesday, December 24, 2014 at 09:45 — #Econometrics