Abstract
The facility location problem (FLP) consists of finding the best location for a facility (service center) to efficiently meet the demands of client points. Studies on the FLP have considered flat or spherical surfaces for computation of distances between these entities. However, the Earth’s surface is not flat or perfectly spherical. This article considers the ellipsoid as a more appropriate model of the Earth’s surface and adapts the FLP to address this assumption. The optimization of the adapted FLP was performed using Solver. Results on geographic instances provided important insights about the logistic financial aspect of the ellipsoid assumption.
Introduction
Strategic planning of the supply chain (SPSC) provides support for the decision-making process related to optimizing a production network and the configuration of the supply chain. 1 Among the decisions involved in SPSC, the following can be mentioned: vertical integration policies, capacity sizing, selection of technology, sourcing, and facility location. 1 Within this context, the facility location problem (FLP) addresses the identification of the optimal supply chain configuration in order to obtain an efficient trade-off between costs, quality, and the availability of products. 1
The general FLP consists of finding the most suitable geographical location to establish (build) service centers or facilities. The decisions regarding the optimal locations involve multiple sites and issues of costs, response times, proximity to certain locations, and accessibility to these sites. 2 Determining the location of the facility is a tactical and long-term investment. 3 As discussed in the study by Shih, 4 “poor facility location decisions can lead to high transportation costs, inadequate supplies of raw materials and labor, loss of competitive advantage, and financial loss. 5 ”
Particularly, decisions about the location of facilities are related to the following issues: The facility must be close to the customers. Priorities are given to keep competitive transportation times. The facility must be located near suppliers of raw materials to reduce costs of materials and labor.
According to ReVelle and Eiselt 6 , the location model has four basic features as follows:
The customers that are located at points or routes.
The facilities that must be located.
The geographic space or surface where the customers and facilities are located.
The metric that indicates the distance, cost, or time between the locations of clients and facilities.
As presented, space is one of the main features for the location model of clients and facilities, and the accuracy of the distances for the facility location model depends on this feature. Most works have considered the space for the FLP to be flat or spherical. 4,7 –10 As discussed in the study by Shih 4 , if the location points are widely separated, the difference in driving distance between the flat surface assumption (i.e. Euclidean distance) and the spherical surface assumption (i.e. great circle distance) may be significant, leading to important variations in the optimal location of the facility. However, as presented in Figure 1, the shape of the Earth, which is the reference space for the FLP, is not flat or perfectly spherical.

Geoid visualization using units of gravity (picture taken from the study by National Aeronautics and Space Administration 11 ). Each point of the geoid’s surface is normal to the direction of gravity.
The geoid, which is considered to represent the truer shape of the Earth, 12 can be approximated as a reference ellipsoid 13 as presented in Figure 2.

Ellipsoid approximation to the geoid.
The ellipsoid of revolution (sometimes called spheroid) has been used as the preferred surface on which point coordinates such as latitude (φ), longitude (λ), and elevation (h) are defined. 14,15
In this work, the surface of the ellipsoid of revolution is analyzed to enhance the mathematical model of the single FLP. This article is structured as follows. In “The ellipsoid and the spherical Earth” section, technical background on the ellipse and the ellipsoid are presented. Then, in “The FLP on ellipsoids” section, the integrated mathematical model of the Weber problem with ellipsoids is presented and discussed. The results of the ellipsoidal model are discussed in “Results” section. Finally, the conclusions and future work are presented in the last section.
The ellipsoid and the spherical Earth
An ellipse is defined as the set of all points P in a plane such that the sum of the distances of P from two fixed points (F and F′) in the plane is constant. 16 Figure 3 presents an overview of the main elements that define an ellipse. Each of the fixed points is called a focus, and together they are called foci. The line segment V′ V through the foci is defined as the major axis, while the perpendicular bisector B′ B of the major axis is defined as the minor axis. Each end of the major axis (i.e. V and V′) is called a vertex. Finally, the midpoint of the line segment F′ F is called the center of the ellipse. 16

Ellipse with center at the origin.
An ellipsoid is the three-dimensional analogue of an ellipse. As presented in Figure 4, the ellipsoid consists of distinct semi-axis lengths a, b, and c. The Earth is approximated to an oblate ellipsoid where a = b > c.

Triaxial ellipsoid with distinct semi-axis lengths.
The standard equation of an ellipsoid centered at the origin of a Cartesian coordinate system is expressed as follows:
Note that the ellipsoid may approximate to a sphere if a = b = c. By knowing this fact, coordinates in the ellipsoid can be expressed in terms of coordinates on the sphere. As presented in Figure 5, the position of a point P on the spherical model of the Earth’s surface is determined by the three main parameters: r (radial distance), θ (longitude), and φ (latitude). 17

Spherical coordinate system.
With these parameters, the parametrization of a sphere of radius r centered at the origin can be expressed as follows 17 –19 :
The ellipsoid can be parametrized by slightly modifying the parametrization of the sphere 19 :
where
As presented in equation (3), the angles θ and φ do not represent geometric angles in the ellipsoid itself; these represent angles in the sphere with r = 1.
Distances on the sphere
The shape of the Earth’s surface is complex, and it is not comparable to any surface that can be modeled with simple mathematical formulation. Under the assumption that the shape of the Earth’s surface can be approximated to a sphere, any i-th point Pi on the surface can be defined by its latitude (φi) and longitude (θi).
The shortest distance or arc length between two points on the spherical surface is defined by:
(a) αij: the angle between the vectors that define the points Pi and Pj from the center of the sphere. 10,20 This angle is estimated as:
(b) R: the radius of the sphere.
The arc length between two points Pi and Pj is expressed as:
The FLP on the sphere
Within the context of the single FLP on the sphere, there are ri(φi, θi) demand points (clients), where i = 1, …, n, and rj(φj, θj) potential locations on which it is possible to place a service point (facility). Each ri(φi, θi) and rj(φj, θj) points are connected by the arcs of length lij, which lie on the spherical surface. A weight wi is associated to each demand point ri(φi, θi). This weight may be associated to important data regarding the demand point such as the magnitude of the demand. 21 Then, the single FLP consists of minimizing the weighted distances from the potential location of the facility to the customers. 8,10,21,22 The objective function is expressed as follows:
The FLP on ellipsoids
Geodesics on the ellipsoid
A geodesic is defined as the shortest path between two points on a curved surface, particularly if the surface is treated as an ellipsoid of revolution.
23
When the Earth is considered as a sphere, the geodesics are great circles (all of which are closed). If the Earth is considered to be an oblate ellipsoid, the equator and the meridians are the only closed geodesics.
24
Figure 6 presents the geodesic AB on an ellipsoid of revolution where
23
: NA and NB are meridians;
s12 is the length of the geodesic AB; λ12 is the longitude of B relative to A;
φ1 and φ2 are the latitudes of A and B, respectively;
α0, α1, and α2 are the azimuths of the geodesic at E, A, and B, respectively.

Geodesic AB on an ellipsoid of revolution where N is the north pole and E, H, and F lie on the equator.
From Figure 6, the arc delimited by the geodesic AB is identified. The geodesic of length s12 connects two points: A at latitude and longitude (φ1 and λ1) and B at (φ2 and λ2). Within this context, two geodesic problems are usually considered 23 :
Direct problem: It consists of finding the end point of a geodesic given its starting point (A), initial azimuth (α1), and length (s12). Thus, B and α2 have to be determined.
Inverse problem: It consists of finding the shortest path between two given points. In this case, s12, α1, and α2 have to be determined.
These problems involve solving the triangle NAB given an angle (α1 for the direct problem and λ12 = λ1 − λ2 for the inverse problem) and two adjacent sides.
For the purposes of the present work, the inverse problem is of particular interest. In the literature, there are methods to calculate the distance between two points on the ellipsoid’s surface as the shortest tangent method. 25 However, the direct and inverse iterative methods of Vincenty 26 are computationally efficient. Also, these methods are appropriate because the Earth’s surface is assumed to be an oblate spheroid, which is a closest approach to the true shape of the Earth. This assumption has led to the use of Vincenty’s methods for partitioning of the ellipsoidal Earth on regions for studies on meteorology, forestry, and biology. 27 In the “ Solution of the inverse problem” section, the basis of the method for the inverse problem of Vincenty 26 is reviewed to describe the adaptation performed for the present work.
Solution of the inverse problem
In “The ellipsoid and the spherical Earth” section, the parametrization of the ellipsoid from spherical coordinates was presented. Vincenty considered this information for the development of the inverse iterative method, and the following notation was defined as 26 :
a and b: the major and minor semi-axes of the ellipsoid (a = 6378137 m according to WGS-84 and b = (1 − f)a where f = 1/298.2572221);
φ: the geodetic latitude, positive north of the equator;
f: the flattening of the spheroid = (a − b)/a;
U: the reduced latitude defined by tanU = (1 − f)tanφ;
λ: the difference in longitude on an auxiliary sphere;
σ: the angular distance between two points, P1 and P2, on the sphere;
S: the ellipsoidal distance between the points P1 and P2.
The pseudocode to estimate the length of s12 (i.e. S) using the inverse method of Vincenty is presented in Figure 7.

Pseudocode to estimate the length of the ellipsoidal distance between two points. 26
Integrated Weber model on ellipsoids
By considering the notation presented in “Solution of the inverse problem” section, the single FLP described by equation (6) can be adapted to the ellipsoid by defining the following terms:
wi: weight associated to each i-th demand point;
sfi: ellipsoidal distance between the service point and the i-th demand point.
Then, the integrated Weber problem consists of finding the location of the facility
which is subject to:
where
Results
Different sets of client points were considered to solve the FLP. For comparison purposes, the spherical model of the FLP was considered to analyze the feasibility of the FLP on ellipsoids. The results are discussed in the following sections.
Test set A: 10 randomly selected client points
For initial testing of the FLP model on ellipsoids, the client points presented in Table 1 were considered. The pseudocode presented in Figure 7 was implemented in Microsoft Excel 2010 for each client point. Then, the tool Solver Excel Solver from Frontline Solvers (R) was used to solve the objective function defined by equation (7) with the restrictions given by equation (8) considering all client points. For comparison purposes, the FLP was solved considering the assumption of the spherical Earth given by equation (6). For all client points, the weights were considered constant and equal to 1.0 (i.e. wi = 1 for i = 1,…,10). Table 2 presents the optimal facility location points for the FLP on the ellipsoid and on the sphere. The value of the objective function represents the sum of all distances (in meters) from the facility to all client points. The location of each client point and the facilities (under spherical and ellipsoid models) are presented in Figure 8.
Latitude (ϕ) and longitude (λ) coordinates for the test set A.
Facility location for the test set A: location under the spherical and ellipsoidal models.

Facility location for the test set A: map locations. (Map data ©2016 Google).
As presented in Table 2, there is a difference of 69,186.72 m (approximately 69 km) between the total distances of the spherical and ellipsoid models. In Table 3, the individual distances from all client points to the located facility are presented. It can be observed that for some points, the distances are longer under the spherical model than for the ellipsoid model. However, for other client points, the distances are longer under the ellipsoid model. By considering the sum of all absolute differences, a total difference of approximately 142 km is estimated between the ellipsoid and the spherical distance models.
Facility location for the test set A: individual spherical and ellipsoidal distances.
The distance differences presented in Table 3 must be considered during budget and delivery time planning due to the significance of some of them (i.e. 22.835 km for Cusco, Peru, 27.338 km for Lusaka, Zambia, and 28.330 km for Hobart, Tasmania). Because most of the distance differences are positive (considering the spherical model as reference), it can be assumed that distances on the spherical model are longer than the distances on the ellipsoid model. However, this assumption changes for specific regions.
Due to this finding, budget and delivery planning for the FLP should consider the ellipsoid model because the spherical model does not consider changes in the Earth’s surface across all regions. The ellipsoid model considers these changes for each Earth’s quadrant. The information presented in Table 3 corroborates this finding.
Test set B: Three client points on the Western Hemisphere and seven on the Eastern Hemisphere
The second testing considered the robustness of the FLP under the spherical and ellipsoid models for changes in the density of client points. For this case, the points presented in Table 4 were considered. Table 5 presents the optimal facility location points for the FLP on the ellipsoid and on the sphere for these points. Then, the location of each client point and the facilities (under spherical and ellipsoid models) are presented in Figure 9.
Latitude (φ) and longitude (λ) coordinates for the test set B.
Facility location for the test set B: location under the spherical and ellipsoidal models.

Facility location for the test set B: map locations (Map data ©2016 Google).
According to the results presented in Table 5 and Figure 9, it can be observed that the facility is to be located within the region where most of the client points are located. This is observed for the FLP under the spherical and ellipsoid models. As presented in Table 5, there is a difference of approximately 17.71 km between the total distance of the ellipsoid model and the spherical model. These results are consistent with those presented in Table 2.
Test set C: Seven client points on the Western Hemisphere and three on the Eastern Hemisphere
The final testing considered an inverse balance on the density of client points presented in “Test set B: Three client points on the Western Hemisphere and seven on the Eastern Hemisphere” section. As presented in Figure 10, the facility (under the spherical and ellipsoid models) is located within the region where most of the client points are located (in this case, in the Western Hemisphere). In Table 6, the location for these facilities is presented where a difference of approximately 43.34 km exists between the total distances of both models.

Facility location for the test set C: map locations (Map data ©2016 Google).
Facility location for the test set C: location under the spherical and ellipsoidal models.
Conclusions and future work
Based on the results presented in “Test set A: 10 randomly selected client points”, “Test set B: Three client points on the Western Hemisphere and seven on the Eastern Hemisphere, and “Test set C: Seven client points on the Western Hemisphere and three on the Eastern Hemisphere” sections, minimum total coverage distance is achieved with the FLP considering the ellipsoid as the model for the Earth’s surface. This distance is smaller than the total coverage distance achieved under the spherical model.
However, depending on the location of the client points across the hemispheres, the magnitude of this difference may increase or decrease. This is because the ellipsoid model considers quadrant-dependent changes in the Earth’s surface. This is an assumption closer to the real surface changes of the Earth which are not considered by the spherical model. Hence, more accurate facility location results can be obtained with the ellipsoid logistic model presented in this work.
These quadrant-dependent distance changes must be considered by the enterprise because budget planning related to factors such as lead times, transportation costs, CO2 emissions, and so on is dependent on the estimated travel distances from the facility to the client points. 28,29 Thus, current cost optimization approaches to strategic logistic costs should incorporate the ellipsoid model for the FLP and routing plans.
Finally, the FLP on the ellipsoid has the same behavior than the FLP on the sphere or on the plain regarding the location of the facility when most of the client points are concentrated within a specific region. In this case, the facility will be located within the region where most of the clients are located. This observation is important to preliminarily estimate the jurisdiction of the facility.
The present work can be extended into the following points regarding the FLP and vehicle routing problem fields:
Footnotes
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship and/or publication of this article: Our institution, Universidad Popular Autonoma del Estado de Puebla, provided the funding for the Article Processing Charge (APC) of 600USD.
