# A Novel Exoplanetary Habitability Score via Particle Swarm Optimization of CES Production Functions

###### Abstract

The search for life has two goals essentially: looking for planets with Earth-like conditions (Earth similarity) and looking for the possibility of life in some form (habitability). Determining habitability from exoplanet data requires that determining parameters are collectively considered before coming up with a conclusion as no single factor alone contributes to it. Our proposed models, would serve as an indicator while looking for new habitable worlds, if computed with precision and efficiency. The models are of the type constrained optimization, multivariate, convex but may suffer from curvature violation and premature convergence impacting desired habitability scores. We mitigate the problem by proposing modified Particle Swarm Optimization (PSO) to tackle constraints and ensuring global optima. In the process, a python library to tackle such problems has been created.

## I Introduction

The search for extra-terrestrial life [1, 2] and potentially habitable extrasolar planets [3, 4] has been an international venture since Frank Drake’s attempt with Project Ozma in the mid-20th century [5]. Cochran, Hatzes, and Hancock[6] confirmed the first exoplanet in 1991. This marked the start of a trend that has lasted 25 years and yielded over 3,700 confirmed exoplanets. There have been attempts to assess the habitability of these planets and to assign a score based on their similarity to Earth. Two such habitability scores are the Cobb-Douglas Habitability (CDH) score [7, 8] and the Constant Elasticity Earth Similarity Approach (CEESA) score. Estimating these scores involves maximizing a production function while observing a set of constraints on the input variables.

Under most paradigms, maximizing a continuous function requires calculating a gradient. This may not always be feasible for non-polynomial functions in high-dimensional search spaces. Further, subjecting the input variables to constraints, as needed by CDH and CEESA, are not always straightforward to represent within the model. This paper details an approach to Constrained Optimization (CO) using the swarm intelligence metaheuristic. Particle Swarm Optimization (PSO) is a method for optimizing a continuous function that does away with the need for calculating the gradient. It employs a large number of randomly initialized particles that traverse the search space, eventually converging at a global best solution encountered by at least one particle [9, 10].

Particle Swarm Optimization is a distributed method that requires simple mathematical operators and short segments of code, making it a lucrative solution where computational resources are at a premium. Its implementation is highly parallelizable. It scales with the dimensionality of the search space. The standard PSO algorithm does not deal with constraints but, through variations in initializing and updating particles, constraints are straightforward to represent and adhere to, as seen in Section III-B. Poli[11, 12] carried out extensive surveys on the applications of PSO, reporting uses in Communication Networks, Machine Learning, Design, Combinatorial Optimization and Modeling, among others.

This paper demonstrates the applicability of Particle Swarm Optimization in estimating habitability scores, CDHS and CEESA of an exoplanet by maximizing their respective production functions (discussed in Sections II-A and II-B). CDHS considers the planet’s Radius, Mass, Escape Velocity and Surface Temperature, while CEESA includes a fifth parameter, the Orbital Eccentricity of the planet. The Exoplanet Catalog hosted by the Planetary Habitability Laboratory, UPR Arecibo records these parameters for each exoplanet in Earth Units [13]. Section V reports the performance of PSO and discusses the distribution of the habitability scores of the exoplanets.

## Ii Habitability Scores

### Ii-a Cobb-Douglas Habitability Score

Estimating the Cobb-Douglas Habitability (CDH) score [7] requires estimating an interior score () and a surface score () by maximizing the following production functions,

R^α. | D^β , | (1a) | |||||

V_e^γ. | T_s^δ , | (1b) |

where, , , and are density, radius, escape velocity and surface temperature respectively. , , and are the elasticity coefficients subject to . Equations 1a and 1b are convex under either Constant Returns to Scale (CRS), when and , or Decreasing Returns to Scale (DRS), when and . The final CDH score is the convex combination of the interior and surface scores as given by,

(2) |

### Ii-B Constant Elasticity Earth Similarity Approach Score

The Constant Elasticity Earth Similarity Approach (CEESA) uses the following production function to estimate the habitability score of an exoplanet,

(3) |

where, is the fifth parameter denoting Orbital Eccentricity. The value of lies within . The coefficients , , , and lie in and sum to , . The value of is constrained by the scale of production used, under DRS and under CRS.

## Iii Particle Swarm Optimization

Particle Swarm Optimization (PSO) [9] is a biologically inspired metaheuristic for finding the global minima of a function. Traditionally designed for unconstrained inputs, it works by iteratively converging a population of randomly initialized solutions, called particles, toward a globally optimal solution. Each particle in the population keeps track of its current position and the best solution it has encountered, called . Each particle also has an associated velocity used to traverse the search space. The swarm keeps track of the overall best solution, called . Each iteration of the swarm updates the velocity of the particle towards its and the values.

### Iii-a PSO for Unconstrained Optimization

Let be the function to be minimized, where is a -dimensional vector. is also called the fitness function. Algorithm 1 outlines the approach to minimizing using PSO. A set of particles are randomly initialized with a position and a velocity, where and are the lower and upper boundaries of the search space. The position of the particle corresponds to its associated solution. The algorithm initializes each particle’s to its initial position. The position that corresponds to the minimum fitness is selected to be the position of the swarm. Shi and Eberhart[10] discussed the use of inertial weights to regulate velocity to balance the global and local search. Upper and lower bounds limit velocity within .

On each iteration, the algorithm updates the velocity and position of each particle. For each particle, it picks two random numbers from a uniform distribution, and updates the particle velocity. Here, is the inertial weight and are the global and particle learning rates. If the new position of the particle corresponds to a better fit than its , the algorithm updates to the new position. Once the algorithm has updated all particles, it updates to the new overall best position. A suitable termination criteria for the swarm, under convex optimization, is to terminate when has not changed by the end of an iteration.

then |

### Iii-B PSO with Leaders for Constrained Optimization

Although PSO has eliminated the need to estimate the gradient of a function, as seen in Section III-A, it still is not suitable for constrained optimization. The standard PSO algorithm does not ensure that the initial solutions are feasible, and neither does it guarantee that the individual solutions will converge to a feasible global solution. Solving the initialization problem is straightforward, resample each random solution from the uniform distribution until every initial solution is feasible. To solve the convergence problem each particle uses another particle’s value, called , instead of its own to update its velocity. Algorithm 2 describes this process.

On each iteration, for each particle, the algorithm first picks two random numbers . It then selects a value from all particles in the swarm that is closest to the position of the particle being updated as its . The value substitutes in the velocity update equation. While updating for the particle, the algorithm checks if the current fit is better than , and performs the update if the current position satisfies all constraints. The algorithm updates as before.

and | |

satisfies all constraints | then |

## Iv Representing the Problem

A Constrained Optimization problem can represented as,

subject to | ||||

Ray and Liew[14] describe a way to represent non-strict inequality constraints when optimizing using a particle swarm. Strict inequalities and equality constraints need to be converted to non-strict inequalities before being represented in the problem. Introducing an error threshold converts strict inequalities of the form to non-strict inequalities of the form . A tolerance is used to transform equality constraints to a pair of inequalities,

l = 1…r , | |||||

l = 1…r . |

Thus, equality constraints become inequality constraints, raising the total number of constraints to . For each solution , denotes the constraint vector where, . When , the solution lies within the feasible region. When , the solution violates the ^{th} constraint.

### Iv-a Representing CDH Score Estimation

### Iv-B Representing CEESA

The representation of CEESA score estimation (described in Section II-B) under DRS is,

d.D^ρ+ t.T_s^ρ | (7) | |||||

v.V_e^ρ+ e.E^ρ) ^ηρ | ||||||

subject to | (8a) | |||||

(8b) | ||||||

(8c) | ||||||

(8d) | ||||||

(8e) | ||||||

(8f) |

Under CRS there is no need for the parameter (since ). Thus, the objective function for the problem reduces to,

d.D^ρ+ t.T_s^ρ | (9) | |||||

v.V_e^ρ+ e.E^ρ) ^1ρ |

## V Experiment and Results

Parameter | Description | Unit |
---|---|---|

P. Radius | Estimated radius | Earth Units (EU) |

P. Density | Density | Earth Units (EU) |

P. Esc Vel | Escape velocity | Earth Units (EU) |

P. Ts Mean | Mean Surface temperature | Kelvin (K) |

P. Eccentricity | Orbital eccentricity |

The data set used for estimating the Habitability Scores of exoplanets was the Confirmed Exoplanets Catalog maintained by the Planetary Habitability Laboratory (PHL) [13]. The catalog records observed and modeled parameters for exoplanets confirmed by the Extrasolar Planets Encyclopedia. Table I describes the parameters from the PHL Exoplanets’ Catalog (PHL-EC) used for the experiment. Since surface temperature and eccentricity are not recorded in Earth Units, we normalized these values by dividing them with Earth’s surface temperature () and eccentricity (). PHL-EC assumes an Eccentricity of 0 when unavailable. The PHL-EC records empty values for planets whose surface temperature is not known. We chose to drop these records from the experiment.

The implementation resulted in a Python library now available on the Python Packaging Index under the name PSOPy [15]. It can be installed through pip as pip install psopy (also available in the github folder AstrIRG). Our implementation used particles to traverse the search space, with learning rates and . It used an inertial weight of and upper and lower bounds . We used an error threshold of to convert strict inequalities to non-strict inequalities, and a tolerance of to transform an equality constraint to a pair of inequalities. Further implementation details are discussed in the Appendix.

The plots in Figures (a)a and (b)b describe the distribution of the CDH scores across exoplanets tested from the PHL-EC. Figures (c)c and (d)d show the distribution of iterations required to converge to a global maxima. The spike at 0 is caused by particles converging to a that does not shift from the original position (for a more detailed explanation see Appendix A). The plots in Figures (b)b and (b)b describe the distribution of the CEESA score across the exoplanets, while Figures (d)d and (c)c show the distribution of iterations to convergence. These graphs aggregate the results of optimizing the Habitability Production Functions (Equations 4, 6a, 7 and 9) for each exoplanet in the PHL-EC by the method described in Algorithm 2.

Table II records the CDH scores for a sample of exoplanets under CRS at and . , , and record the parameters of Equation 4. and record the maxima for the objective functions. and specify the number of iterations taken to converge to a stable value. Under the Class column there are four categories for the planets — Psychroplanets (psy), Mesoplanets (mes), Non-Habitable planets (non) and Hypopsychroplanets (hyp). Table III records the CDH scores for a sample under DRS, with , , and recording the parameters of Equation 4.

Tables IV and V record the estimated CEESA scores under CRS and DRS respectively. , , , , , and record the parameters of Equation 7 in Table V and the parameters of Equation 9 in Table IV. However, since under the CRS constraint, , there is no need for the parameter in Table IV. specifies the number of iterations taken to converge to the maxima.

These tables indicate that although CEESA has 7 parameters and 16 constraints under DRS, PSO takes a little over twice the number of iterations to converge as in each step of the CDH score estimation, which has 2 parameters and 5 constraints. This is a promising result as it indicates that the iterations required for converging increases sub-linearly with the number of parameters in the model. As for real time taken to converge, PSO took () to estimate the CDH score under CRS for exoplanets, at an average of for each planet for each individual score (interior and surface) of the CDH score. For CDH estimation under DRS, it took () at an average of for each part of the CDH score. The CEESA calculations, requiring a single estimate, took a little over half the CDH estimation execution time to run. Under DRS it took a total of () at per planet, while under CRS it took () at per planet.

Name | Class | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

GJ 176 b | non | |||||||||

GJ 667 C b | non | |||||||||

GJ 667 C e | psy | |||||||||

GJ 667 C f | psy | |||||||||

GJ 3634 b | non | |||||||||

HD 20794 c | non | |||||||||

HD 40307 e | non | |||||||||

HD 40307 f | non | |||||||||

HD 40307 g | psy | |||||||||

Kepler-186 f | hyp | |||||||||

Proxima Cen b | psy | |||||||||

TRAPPIST-1 b | non | |||||||||

TRAPPIST-1 c | non | |||||||||

TRAPPIST-1 d | mes | |||||||||

TRAPPIST-1 e | psy | |||||||||

TRAPPIST-1 g | hyp |

Name | Class | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

GJ 176 b | non | |||||||||

GJ 667 C b | non | |||||||||

GJ 667 C e | psy | |||||||||

GJ 667 C f | psy | |||||||||

GJ 3634 b | non | |||||||||

HD 20794 c | non | |||||||||

HD 40307 e | non | |||||||||

HD 40307 f | non | |||||||||

HD 40307 g | psy | |||||||||

Kepler-186 f | hyp | |||||||||

Proxima Cen b | psy | |||||||||

TRAPPIST-1 b | non | |||||||||

TRAPPIST-1 c | non | |||||||||

TRAPPIST-1 d | mes | |||||||||

TRAPPIST-1 e | psy | |||||||||

TRAPPIST-1 g | hyp |

Name | Class | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

GJ 176 b | non | |||||||||

GJ 667 C b | non | |||||||||

GJ 667 C e | psy | |||||||||

GJ 667 C f | psy | |||||||||

GJ 3634 b | non | |||||||||

HD 20794 c | non | |||||||||

HD 40307 e | non | |||||||||

HD 40307 f | non | |||||||||

HD 40307 g | psy | |||||||||

Kepler-186 f | hyp | |||||||||

Proxima Cen b | psy | |||||||||

TRAPPIST-1 b | non | |||||||||

TRAPPIST-1 c | non | |||||||||

TRAPPIST-1 d | mes | |||||||||

TRAPPIST-1 e | psy | |||||||||

TRAPPIST-1 g | hyp |

Name | Class | |||||||||
---|---|---|---|---|---|---|---|---|---|---|

GJ 176 b | non | |||||||||

GJ 667 C b | non | |||||||||

GJ 667 C e | psy | |||||||||

GJ 667 C f | psy | |||||||||

GJ 3634 b | non | |||||||||

HD 20794 c | non | |||||||||

HD 40307 e | non | |||||||||

HD 40307 f | non | |||||||||

HD 40307 g | psy | |||||||||

Kepler-186 f | hyp | |||||||||

Proxima Cen b | psy | |||||||||

TRAPPIST-1 b | non | |||||||||

TRAPPIST-1 c | non | |||||||||

TRAPPIST-1 d | mes | |||||||||

TRAPPIST-1 e | psy | |||||||||

TRAPPIST-1 g | hyp |

## Vi Conclusion

Particle Swarm Optimization mainly draws its advantages from being easy to implement and highly parallelizable. The algorithms described in Section III use simple operators and straightforward logic. What is especially noticeable is the lack of the need for a gradient, allowing PSO to work in high dimensional search spaces with a large number of constraints, precisely what is needed in a potential Habitability score estimate. Further, particles of the swarm, in most implementations operate independently during each iteration, their updates can occur simultaneously and even asynchronously, yielding much faster execution times than those outlined in Section V. However, since strict inequalities and equality constraints are not exactly represented, the resulting solution may not be as accurate as direct methods. Despite this, using PSO to calculate the habitability scores is beneficial when the number of input parameters are large, which further increases the number of constraints, resulting in a model too infeasible for traditional optimization methods.

Determining habitability from exoplanet requires that determining parameters are collectively considered [13] before coming up with a conclusion as no single factor alone contributes to it. Our proposed model would serve as an indicator while looking for new habitable worlds. Eccentricity may have some effect on habitabilty and the models for computation should address that. CDHS doesn’t, at least for the Trappist system (otherwise considered a set of potentially habitable exoplanets) since the eccentricities for all memebers of the Trappist system are 0 identically. CDHS, being a product model then will render the habitability score to be 0, not in agreement with observations and overall opinion in the community. This is the reason CESSA is considered in the process of habitability score computation (additive nature of the model). One might wonder why metaheuristic optimization was applied on two different optimization problems. We hope, our clarification would suffice.

However, the functional forms considered to compute the habitability score pose challenges. As we intend to add more parameters (such as eccentricity) to the basic model [13][7], the functional form tends to suffer from curvature violation [16][17]. Even though global optima is guaranteed, premature convergence and local oscillations are hard to mitigate. An attempt to address such issues, with moderate success, could be found in [18]. The greatest contribution of the manuscript is to propose an evolutionary algorithm to track dynamic functions of the type that allow for the oscillation that were instead mitigated with SGA in [17]. Consequently, a Pyhton library is integrated with the open source tool suite, an add on for coding enthusiasts to test our method.

## Appendix A Ensuring Convergence

Although the termination criteria illustrated in Algorithm 1, , is a lucrative solution for unconstrained optimization, it might cause early termination for a constrained problem. This could occur in Algorithm 2 as the particles traverse the search space. Due to the way they are updated, all individual particles could move to infeasible regions. There are no updates made to their corresponding s, leading to no change in . Thus, although the swarm has not yet converged to the global optima, it terminates and reports the current as the optima.

To prevent this from happening our implementation maintained a counter that it incremented every time changed by a value less than the threshold. When the counter reached , the algorithm terminated and reported the current value of as the global optima. If the value of changed by more than the threshold, it reset the counter to . Since the function is convex, and all points are initialized to feasible solutions (thus, ensuring that a feasible has been encountered) there was no need to restrict the total number of possible iterations.

This explains the spike at in Figure (c)c. For every exoplanet, the implementation waits a iterations for declaring convergence, thus overstating the total iterations by . To adjust for this offset, we subtracted from the total iterations before constructing the histogram plots in Figures (c)c, (d)d, (c)c and (d)d. Therefore, although the graphs report iterations to convergence for a significant number of planets, only has not shifted, the particles have actually converged to a common point around a region within the threshold of the global best.

## Appendix B Improving Performance

Most programming languages today have either built-in or external library support for linear algebra routines. The implementation of these libraries are well-tested and efficient, and take advantage of system features such as multiple cores or a general purpose GPU computing device. Taking advantage of these libraries requires representing most entities as either matrices or vectors.

Matrices and represent current position and velocity of all particles in the swarm.

∣ | ∀i=1…n, ∀j=1…d ] , | |||||

∣ | ∀i=1…n, ∀j=1…d ] . |

and denote position and velocity of particle in dimension . Matrices and represent the and values for the swarm, where,

Section IV outlined the use of a constraint vector for describing the feasibility of a solution. This concept generalizes to a matrix , where,

Let be two random vectors of length drawn from the uniform distribution . We use the shorthand notation of to denote the ^{th} row of some matrix . Let be the fitness function and construct for solutions . denotes the matrix and denotes the Hadamard product. The update of each particle in Algorithm 2 is implemented as,

Finally, after all particles are updated, is updated according to,

## Acknowledgment

The authors would like the Science and Engineering Research Board (SERB)-Department of Science and Technology (DST), Government of of India for supporting this research. The project reference number is: SERB-EMR/ 2016/005687.

## References

- [1] S. Shostak. (2017, Dec.) Do aliens exist? your tax dollars may hold the truth. [Online]. Available: http://www.newsweek.com/do-aliens-exist-your-tax-dollars-may-hold-answer-753151
- [2] E. W. Schwieterman et al., “Exoplanet biosignatures: A review of remotely detectable signs of life,” arXiv preprint arXiv:1705.05791, May 2017.
- [3] M. Johnson. (2018, Jan.) Kepler and K2 mission overview. [Online]. Available: https://www.nasa.gov/mission_pages/kepler/overview
- [4] M. C. LoPresto and H. Ochoa, “Searching for potentially habitable extra solar planets: a directed-study using real data from the NASA kepler-mission,” Physics Education, vol. 52, no. 6, p. 065016, 2017.
- [5] H. P. Shuch, Searching for Extraterrestrial Intelligence: SETI past, present, and future. Springer Science & Business Media, 2011, ch. 2, pp. 13–18.
- [6] W. D. Cochran, A. P. Hatzes, and T. J. Hancock, “Constraints on the companion object to hd 114762,” The Astrophysical Journal, vol. 380, pp. L35–L38, Oct. 1991.
- [7] K. Bora, S. Saha, S. Agrawal, M. Safonova, S. Routh, and A. Narasimhamurthy, “Cd-hpf: New habitability score via data analytic modeling,” Astronomy and Computing, vol. 17, pp. 129–143, 2016.
- [8] S. Saha, S. Basak, K. Bora, M. Safonova, S. Agrawal, P. Sarkar, and J. Murthy, “Theoretical validation of potential habitability via analytical and boosted tree methods: An optimistic study on recently discovered exoplanets,” arXiv preprint arXiv:1712.01040, Dec. 2017.
- [9] R. Eberhart and J. Kennedy, “A new optimizer using particle swarm theory,” in Micro Machine and Human Science, 1995. MHS’95., Proceedings of the Sixth International Symposium on. IEEE, 1995, pp. 39–43.
- [10] Y. Shi and R. Eberhart, “A modified particle swarm optimizer,” in Evolutionary Computation Proceedings, 1998. IEEE World Congress on Computational Intelligence., The 1998 IEEE International Conference on. IEEE, 1998, pp. 69–73.
- [11] R. Poli, “An analysis of publications on particle swarm optimization applications,” Essex, UK: Department of Computer Science, University of Essex, 2007.
- [12] ——, “Analysis of the publications on the applications of particle swarm optimisation,” Journal of Artificial Evolution and Applications, 2008.
- [13] A. MÃ©ndez. (2017, Nov.) Phl’s exoplanets catalog. [Online]. Available: http://phl.upr.edu/projects/habitable-exoplanets-catalog/data/database
- [14] T. Ray and K. M. Liew, “A swarm with an effective information sharing mechanism for unconstrained and constrained single objective optimisation problems,” in Evolutionary Computation, 2001. Proceedings of the 2001 Congress on, vol. 1. IEEE, 2001, pp. 75–80.
- [15] A. J. Theophilus, S. Saha, and S. Basak, “PSOPy: A Python implementation of Particle Swarm Optimization,” 2018–, [Online; accessed May 23, 2018]. [Online]. Available: https://pypi.org/project/psopy/
- [16] S. Saha, J. Sarkar, A. Dwivedi, N. Dwivedi, A. M. Narasimhamurthy, and R. Roy, “A novel revenue optimization model to address the operation and maintenance cost of a data center,” Journal of Cloud Computing, vol. 5, no. 1, p. 1, 2016.
- [17] J. Sarkar, B. Goswami, S. Saha, and S. Kar, “Cdsfa stochastic frontier analysis approach to revenue modeling in large cloud data centers,” arXiv preprint arXiv:1610.00624, 2016.
- [18] S. Agrawal, S. Basak, S. Saha, K. Bora, and J. Murthy, “A comparative analysis of the cobb-douglas habitability score (cdhs) with the earth similarity index (esi),” arXiv preprint arXiv:1804.11176, 2018.