以遺傳基因演算法為基礎的符號式回歸引擎及其於全球定位系統之座標轉換上的應用

全文

(1)國立高雄大學電機工程研究所碩士論文. 以遺傳基因演算法為基礎的符號式回歸引擎及其於全球定位系統之座標轉換上的應用 A GP-Based Symbolic Regression Engine for the Transformation of GPS Coordinates. 研究生：周鴻儒撰指導教授：吳志宏博士. 中華民國九十六年七月.

(2)

(3) 致謝本論文能夠完成要感謝許多人，第一位當然非指導教授莫屬，吳志宏老師幫助我完成了人生的第一本著作，無論是在尋找問題、研究方法、實驗架構、程式設計和文法寫作等各方面，吳老師都給予我非常嚴謹的指導，受益良多，除了專業技能之外，透過三年來的相處，也學習了待人處世的觀念和技巧，真的很感謝吳老師。實驗室的學長姐，千慧學姊、吉原學長、瓊輝學長、獻仁學長、唐龍學長、唐虎學長，從入學到現在，受到很多照顧，讓我可以專心順利的在研究工作上。同儕門，香君、宇傑，前前後後我們都畢業了，這段時間我們彼此扶助，互相鼓勵，成為了我繼續學術這條路的動力。碩士班的學弟妹，韋瀚、瀞誼，感謝你們在行政事務、實驗室事務與我ㄧ起分擔責任和工作。大學部的學弟，志恩、銘宏、佳霖，以及其他 e-Eureka 志工團隊的夥伴，身為學長常常請你們幫忙，謝謝你們。學佳，與你在一起兩年了，謝謝妳陪我分擔研究的酸甜苦辣，妳也幫了我好多，與你在一起是吳老師自稱最偉大的成就，希望未來能長長久久走下去。感謝洪宗貝老師、林文揚老師、吳俊興老師、賴智錦老師，除了課業的教授外，在其他的時間，也願意提供學生建議，對於學生能夠考取博士班，順利完成碩士學位，各位老師的幫忙，我非常感謝。我還要感謝學校裡其他的高大人，電機系系辦助理，姿蓉、系主任，徐忠枝老師、資工系系辦助理，淑真、教務處的大哥大姐，總務處和會計室的大哥大姐，過去一段時間的幫助，我銘記在心。感謝父母親和家人，能夠支持我繼續學術這條充滿考驗的路，大家過去的幫助以及祝福，讓我取得碩士學位，未來我相信我可以順利成為鴻儒博士。這篇論文不是我完成的，而是大家！謝謝大家！. 鴻儒. 誌於. 國立高雄大學理工第一實驗大樓智慧計算與應用實驗室民國 96 年 7 月 31 日.

(4) 以遺傳基因演算法為基礎的符號式回歸引擎及其於全球定位系統之座標轉換上的應用指導教授：吳志宏博士國立高雄大學電機工程研究所. 學生：周鴻儒國立高雄大學電機工程研究所. 摘要. 全球定位系統之三維座標需要透過複雜的階層式數學計算才可轉換為二維座標。但是微小的誤差不可避免的會在計算中被累積，使得精準度降低。在本論文中，我們提出一種新的符號式回歸引擎來產生 GPS 座標轉換的回歸方程式。回歸方程式是透過樹狀結構所定義，再藉由遺傳基因演算法的交配、突變與選擇來改變回歸方程式的組合來尋找合適的方程式組合。為了增加回歸引擎的效率，我們採用了可動態改變機率的機制，加上樹狀修飾來降低公式複雜度，但是仍維持原公式的功能。若能夠得到計算複雜度低的座標轉換公式，就能夠使得定位裝置在獨立電源的環境中，延長待機時間。此方法雖然不能夠降低 GPS 的系統誤差，但是透過定位資料成功地降低此 GPS 的統計誤差。此研究的實驗資料採用 TWD67 的大地座標系統，實驗結果顯示，我們確實能夠使用自行設計的符號式回歸引擎來取得計算量低的 GPS 座標轉換公式。關鍵字：遺傳基因演算法、符號式回歸、全球定位衛星系統、座標轉換、TWD67.

(5) A GP-Based Symbolic Regression Engine for the Transformation of GPS Coordinates Advisor(s): Dr. Chih-Hung Wu Department of Electrical Engineering National University of Kaohsiung. Student: Hung-Ju Chou Department of Electrical Engineering National University of Kaohsiung. ABSTRACT Transformation of coordinates usually invokes level-wised processes wherein several sets of complicated equations are calculated. Unfortunately, the accuracy may be corrupted due to the accumulation of inevitable errors between the transformation processes. In this thesis, we present a new symbolic regression engine for generating regressive models for direct transformation from GPS signals to 2-D coordinates. The functional models to be found by the regression engine are represented in tree-based chromosome. Genetic operators like mutation, crossover, and reproduction are applied on such trees. In order to improve the efficiency of the regression process, we design some adaptive mechanisms which adjust the probabilities of applying genetic operations. Methods of trimming and pruning the regressive models for reducing their complexity are also design and implemented. Since target coordinates for a GPS application can be obtained by using simpler transformation formulas, the computational costs and inaccuracy can be reduced. The proposed method, though does not exclude systematic errors due to the imperfection on defining the reference ellipsoid and the reliability of GPS receivers, effectively reduces the statistical errors when the accurate Cartesian coordinates are known from the independent sources. From the experimental results where the target datum TWD67 is investigated, it seems that the proposed method can serve as a direct and feasible solution to the transformation of GPS coordinates. Keywords: GP, symbolic regression, transformation of GPS coordinates, TWD67.

(6) A GP-Based Symbolic Regression Engine for the Transformation of GPS Coordinates by Hung-Ju Chou [email protected] Department of Electrical Engineering, National University of Kaohsiung Kaohsiung City, Taiwan, 811, R.O.C. Advisor: Dr. Chih-Hung Wu [email protected] Department of Electrical Engineering, National University of Kaohsiung Kaohsiung City, Taiwan, 811, R.O.C..

(7) Contents 1 Introduction. 6. 1.1. Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. 1.2. Motivations and Methodologies . . . . . . . . . . . . . . . . . . . . . . .. 7. 1.3. Thesis Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 8. 2 Background and Related Work. 9. 2.1. Transformation of Coordinates . . . . . . . . . . . . . . . . . . . . . . . .. 9. 2.2. Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 14. 3 The Design Concepts. 16. 3.1. Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 16. 3.2. Tree-based Representation . . . . . . . . . . . . . . . . . . . . . . . . . .. 18. 3.3. Defining Fitness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 21. 3.4. Genetic Operations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 3.4.1. Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 22. 3.4.2. Mutation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 24. 3.4.3. Adaptive Selection . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 3.4.4. Tree Refinements . . . . . . . . . . . . . . . . . . . . . . . . . . .. 29. 3.5. Termination Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 31. 3.6. Processing Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 33. 4 Implementation Issues. 35. 4.1. Data Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 35. 4.2. Implementing Genetic Mechanisms . . . . . . . . . . . . . . . . . . . . .. 36. 4.2.1. The Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . .. 38. 4.2.2. The Selection Mechanism . . . . . . . . . . . . . . . . . . . . . .. 38. 4.2.3. The Crossover Mechanism . . . . . . . . . . . . . . . . . . . . . .. 40. 4.2.4. The Mutation Mechanism . . . . . . . . . . . . . . . . . . . . . .. 41. 4.2.5. The Improvement Mechanism . . . . . . . . . . . . . . . . . . . .. 47. 5 Experiments. 58. 5.1. Data Collection and Performance Measurements . . . . . . . . . . . . . .. 58. 5.2. Experiment–A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 60. 5.2.1. 60. Data collection and experiment design . . . . . . . . . . . . . . .. 1.

(8) 5.3 5.4. 5.2.2. Experiment–A–I. . . . . . . . . . . . . . . . . . . . . . . . . . . .. 60. 5.2.3. Experiment–A–II . . . . . . . . . . . . . . . . . . . . . . . . . . .. 62. Experiment–B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 66. 5.3.1. Data collection and experiment design . . . . . . . . . . . . . . .. 66. Experiment–C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 75. 5.4.1. Data collection and experiment design . . . . . . . . . . . . . . .. 75. 5.4.2. Experiment–C–I . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 75. 5.4.3. Experiment–C–II . . . . . . . . . . . . . . . . . . . . . . . . . . .. 75. 5.4.4. Experiment–C–III . . . . . . . . . . . . . . . . . . . . . . . . . . .. 76. 6 Discussion and Conclusion. 85. A Algorithms. 86. B Abbreviation list. 95. 2.

(9) List of Figures 2.1. Level-wised transformation from WGS84 to TM2 in Taiwan . . . . . . .. 13. 2.2. Direct transformation from WGS84 to TM2 . . . . . . . . . . . . . . . .. 14. 3.1. Example of chromosome in tree-based representation, nodes are represented by locus:allele. inv means invalid node . . . . . . . . . . . . . . .. 19. 3.2. Example of crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 23. 3.3. Examples of mutation . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 25. 3.4. Examples of selecting candidate node . . . . . . . . . . . . . . . . . . . .. 28. 3.5. Examples of tree modification . . . . . . . . . . . . . . . . . . . . . . . .. 30. 3.6. Constants in the engine . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 32. 3.7. The proposed genetic-based process of finding regressive model . . . . . .. 34. 4.1. Tree-based chromosome: its data structure: (x + y) × (2 − y) . . . . . . .. 35. 4.2. Data structures of tree pools . . . . . . . . . . . . . . . . . . . . . . . . .. 36. 4.3. Overview of symbolic regression engine . . . . . . . . . . . . . . . . . . .. 39. 4.4. Process steps of crossover. . . . . . . . . . . . . . . . . . . . . . . . . . .. 41. 4.5. Steps of crossover-2p procedure . . . . . . . . . . . . . . . . . . . . . . .. 42. 4.6. CPa is selected by level-based selection, but CPb isn’t . . . . . . . . . . .. 45. 4.7. CPb is selected by level-based selection, but CPa isn’t . . . . . . . . . . .. 45. 4.8. CPa and CPb are selected by level-based selection . . . . . . . . . . . . .. 47. 4.9. Process steps of mutation . . . . . . . . . . . . . . . . . . . . . . . . . .. 48. 4.10 Examples of mutation . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 49. 5.1. The experiment data in experiment–B . . . . . . . . . . . . . . . . . . .. 61. 5.2. The average and maximum deviations between tE and TE . . . . . . . . .. 64. 5.3. The average and maximum deviations between tN and TN. . . . . . . . .. 65. 5.4. The deviations of the regression functions with sampled data sorted . . .. 65. 5.5. The deviations of the regression functions with sampled data unsorted . .. 66. 5.6. Reference points in sampling area 1 . . . . . . . . . . . . . . . . . . . . .. 67. 5.7. Reference points in sampling area 2 . . . . . . . . . . . . . . . . . . . . .. 68. 5.8. Generating data in sampling area 3 . . . . . . . . . . . . . . . . . . . . .. 69. 5.9. Collecting data in sampling area 4 . . . . . . . . . . . . . . . . . . . . . .. 70. 3.

(10) 5.10 Advanced analysis on the result of Experiment–C using different numbers of sampling points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 74. 5.11 The standard reference points partition by geographical distribution . . .. 77. 5.12 The standard reference points partition by height distribution . . . . . .. 78. 5.13 The distribution of sampled data and standard reference points . . . . .. 82. 5.14 Experimental results on synthetic datasets–Results of E1–E4 on TM2E .. 82. 5.15 Experimental results on synthetic datasets–Results of E1–E4 on TM2N .. 83. 5.16 Experimental results on synthetic datasets–Results of N1–N4 on TM2E .. 83. 5.17 Experimental results on synthetic datasets–Results of N1–N4 on TM2N .. 83. 4.

(11) List of Tables 2.1. Parameters used in Eq.(2.1)–Eq.(2.15) : Parameters defined in WGS84 and TWD67 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 11. 2.2. Parameters used in Eq.(2.1)–Eq.(2.15) : More specific parameters . . . .. 12. 3.1. Definition of operators in tree structure . . . . . . . . . . . . . . . . . . .. 20. 3.2. All defined cases for modification . . . . . . . . . . . . . . . . . . . . . .. 29. 4.1. The definitions of structures of a tree . . . . . . . . . . . . . . . . . . . .. 36. 4.2. Parameters defined in tree pools . . . . . . . . . . . . . . . . . . . . . . .. 37. 4.3. The definitions of structures of a tree . . . . . . . . . . . . . . . . . . . .. 38. 5.1. Parameter settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 59. 5.2. Experimental results of tE in B–I and B–II . . . . . . . . . . . . . . . . .. 62. 5.3. Experimental results of tN in B–I and B–II . . . . . . . . . . . . . . . . .. 63. 5.4. Experimental results of B–I . . . . . . . . . . . . . . . . . . . . . . . . .. 71. 5.5. Experimental results of B–II . . . . . . . . . . . . . . . . . . . . . . . . .. 72. 5.6. Experimental results of B–III . . . . . . . . . . . . . . . . . . . . . . . .. 73. 5.7. Experimental results–using all standard reference points for regression . .. 76. 5.8. Experimental results on partitioned data by their geographic locations . .. 79. 5.9. Best functions on partitioned data by their geographic locations . . . . .. 80. 5.10 Experimental results on partitioned data by their height . . . . . . . . .. 80. 5.

(12) Chapter 1 Introduction. 1.1. Problem Definition. The Global Positioning System (GPS) [1, 40] is a satellite-based navigation system made up of a network of 24 satellites which broadcasts precise timing signals by radio to GPS receivers, allowing users to determine their locations on Earth. With the popularity of general purposed GPS receivers as consumer electronics, GPS has been emerging as a convenient tool for positioning and navigation. Positioning by GPS is carried out in three-dimensional geocentric Cartesian coordinates, X, Y , and Z, which is calculated and transmitted to the receivers as coded information. Signals from a general purposed GPS receiver are usually encoded according to the NMEA-0183 standard, carrying positioning information (longitude, latitude, altitude, etc.) in WGS84 [2]. As in practical applications of positioning or navigation, it is frequently required to convert the coordinates derived in one geographic coordinate system to values expressed in another. Coordinates obtained from a GPS receiver also need to be converted if the target coordinates are described in a different way. Transformations between various coordinate systems involve not only complex, usually non-linear, algebraic formulas, but also some very specific numerical parameters which generally have been established as national and continental geodetic datums. If necessary, the transformation process performs level-wise conversions, starting from the original coordinates to the target coordinates, where each level of conversion takes care of a mapping from one coordinate system to another. Unfortunately, some transformation formulas for specific areas are based on estimation and inherently of errors; and therefore, positioning accuracy is degraded as such errors are accumulated during the transformation process. The more steps for transformations, the more inaccuracies. From the aspect of compu6.

(13) tational costs, the more steps for transformations may derive more power consumption; which is a serious problem when GPS is applied on a mobile device.. 1.2. Motivations and Methodologies. This study tries to eliminate the levels needed for the transformation of coordinates and to produce a condensed process for GPS positioning. Conceptually, the level-wise transformation process can be viewed as the composition of all transformation functions. There seemingly exists a set of transformation equations between the signals received from GPS and the target coordinates used in an application. This intuitively forms a problem of regression. If we can derive such formulas that perform fewer levels of transformation, the inaccuracies and computational costs may be reduced. In this paper, we propose a method based on the techniques of soft-computing for deriving simpler formulas for direct transformation of coordinates from WGS84 to specific coordinate systems. Experiments on target datums TWD67 [3, 4, 5] are conducted. The experimental results show that the proposed method can serve as a direct and feasible solution to the transformation of GPS coordinates. Notably, GPS-based positioning needs to consider at least two types of inaccuracy. Firstly, systematic errors exist due to the imperfection on defining the reference ellipsoid [6], though relative positioning with respect to standard reference points defined in the local datums is accurate and acceptable. Secondly, GPS has several sources of errors, such as the radio signal corruption caused by ionospheric delay, tropospheric delay, satellite clock and receiver clock offsets, receiver noise, receiver calibration and multipath, and the synchronization of positioning data from different tracking stations. Approaches that can partially reduce such errors and provide accurate coordinates for positioning have been proposed [7, 8, 9, 10, 11]. Eliminating such errors mentioned above is not the scope of this paper. Moreover, this study is not to find a single transformation method for positioning in all reference systems. The proposed method has been successfully applied, but not restrictedly, to the transformation of TWD67. The same procedure can be applied to the transformation of other target coordinate systems, provided that accurate Cartesian coordinates are known from independent sources. In the following, we present how symbolic regression and genetic programming cope with the problem of coordinate transformations. Genetic algorithms [12] simulate the natural selection rule, the-best-can-fit, and have been successfully applied to many optimization problems. A genetic algorithm iteratively performs operations of crossover and mutation on the expressions of the target problem and evaluates each new solution with respect to the error function for all input data; until the one with lowest error is accepted. Genetic-based methods are usually ap-. 7.

(14) plied for solving complex problems where no effective algorithm or optimal solution exists. Though genetic-based methods provide approximated solutions, they usually do not need to search all solutions exhaustively and can produce satisfactory solutions [13, 14, 15]. To our best knowledge, this study is the first attempt to use genetic-based methods in the coordinate transformation problems for GPS positioning.. 1.3. Thesis Organization. Organization of the thesis is as follows. In Chapter 2, the process of transformation from GPS signals to 2-D coordinates is briefly introduced and giving a brief review on the related works. Then, we formulate the level-wise transformation process as a regression problem and define the related error functions in Chapter 3. The proposed genetic-based approach and associated methods for improving its efficiency are presented in Chapter 3. The implementation is following shown in Chapter 4. The results of experiments on the proposed method are given in Chapter 5. Finally, Chapter 6 gives the discussion and concludes the work.. 8.

(15) Chapter 2 Background and Related Work. 2.1. Transformation of Coordinates. In this paper, the geodetic datum, TWD67 [3, 41], which is based on GRS67, is employed as the targeted 2-D coordinate system. To project 3-D coordinates onto a 2-D local map, 2◦ -Zone Transverse Mercator projection (TM2) is employed in conjunction with TWD67. For positioning on a 2-D map by using GPS, coordinates in WGS84 have to be transformed into the ones in TM2. Most general purposed GPS receivers support geodetic outputs of latitude(φ84 ), longitude(λ84 ), and the height(h84 ) above the reference surface in WGS84. When working with the transformation, geodetic coordinates in WGS84 are firstly transformed into geographic coordinates (X84 , Y84 , Z84 ) using the following standard transformation equations. 0 X84 = (N84 + h84 ) · cos φ84 · cos λ84. (2.1). Y84 = (N 0 + h84 ) · cos φ84 · sin λ84. (2.2). Z84 = (N 0 · (1 − e2 ) + h84 ) · sin φ84. (2.3). where a84 , (1 − e2 sin2 φ84 )1/2 = a84 · (1 − f84 ), (a284 − b284 )1/2 = , a84. 0 N84 =. b84 e84. (2.4) (2.5) (2.6). a84 , b84 , and f84 are the semi-major axis, semi-minor axis, and the flattening ratio of the reference ellipsoid, WGS84, respectively. 9.

(16) Secondly, the values of Cartesian coordinates, X84 , Y84 , and Z84 are transformed according to the local datum, GRS-67. Several transformation methods such as [16, 17, 18] can be used in this stage. Here we use the Molodensky-Badekas’s formulas. . X67. . . δX. . . Xm. . . X84 − Xm. .          Y67  =  δY  +  Ym  + (1 + s) · R ·  Y84 − Ym  Z67 δZ Zm Z84 − Zm In Eq.(2.7),. . 1.  R =  −εZ. +εZ. +εY. 1 −εX. −εY. (2.7). .  +εX  , 1. Xm , Ym , and Zm are the gravity center of WGS84. Then, the obtained GRS67(X67 , Y67 , Z67 ) is converted back to GRS67(φ67 , λ67 , h67 ) by the formulas below. Z67 + ε2 · b67 · sin3 µ67 ) P − e267 · a67 · cos3 µ67 Y67 = tan−1 ( ) X67 P − N0 = cos φ67. φ67 = tan−1 (. (2.8). λ67. (2.9). h67. (2.10). where q P =. 2 X67 + Y672 ,. Z67 · a67 ), P · b67 (a267 − b267 )1/2 ε = , b67. µ = tan−1 (. (2.11) (2.12) (2.13). 0 and b67 , N67 , and e67 are computed using the same formulas in Eq.(2.5), Eq.(2.4), and. Eq.(2.6), respectively, with proper parameters substituted. T M2 Finally, the values of GRS67(φ67 , λ67 , h67 ) is projected to 2-D local eastings E67 T M2 by 2◦ -Zone Transverse Mercator projection (TM2) [19] as follows. and northings N67 0 T M2 · [δλ · cos φ67 + = Ws + K0 · N67 E67 3 3 δλ · cos φ67 · (1 − t2 + η 2 ) + 6 δλ5 · cos5 φ67 · (5 − 18t2 + t4 + 14η 2 − 120 58t2 η 2 + 13η 4 + 4η 6 − 67t2 η 4 − 24t2 η 6 ) + δλ7 · cos7 φ67 · (61 − 479t2 + 179t4 − t6 ) ] 5040. 10. (2.14).

(17) Table 2.1: Parameters used in Eq.(2.1)–Eq.(2.15) : Parameters defined in WGS84 and TWD67 Datum. WGS84. TWD67. Ref. Ellipsoid. WGS84. GRS67. Major Axis. 6378137 (a84 ). 6378160 (a67 ). Minor Axis. 6356752.3142 (b84 ). 6356774.7192 (b67 ). Flattening. 1/298.257223563 (f84 ). 1/298.25 (f67 ). 1st Eccentricity. 0.081819190928906(e84 ). 0.081820179987095 (e67 ). 2nd Eccentricity. 0.082094438036854(ε84 ). 0.082095437110624(ε67 ). · T M2 N67. 0. = K0 · N ·. Sφ δλ2 + · sin φ67 · cos φ67 + N0 2. (2.15). δλ4 · sin φ67 · cos3 φ67 · (5 − t2 + 9η 2 + 4η 4 ) + 24 δλ6 · sin φ67 · cos5 φ67 · (61 − 58t2 + t4 + 720 270η 2 − 330t2 η 2 + 445η 4 + 324η 6 − 680t2 η 4 + 88η 8 − 600t2 η 6 − 192t2 η 8 ) + δλ8 · sin φ67 · cos7 φ67 · (1385 − 311t2 + 40320 543t4 − t6 ) ] 0 In the above formulas, δλ=λ67 − λ0 , t=tan φ67 , and η=ε67 · cos φ67 , where N67 and ε67 are. computed by the formulas in Eq.(2.4) and Eq.(2.13), respectively, with proper parameters substituted. Specific parameters involved in Eq.(2.1)–Eq.(2.15) referring to the reference ellipsoids are given in Table 2.1 and Table 2.2. [41, 42, 43] Level-wise transformation, as shown in Figure 2.1, has been used for a long time and many software tools are developed based on these formulas. Several problems observed from such processes are as follows. • Clearly, Eq.(2.1)–Eq.(2.15) not only include complex, non-linear, high-order, transcendental formulas, but also specific numerical parameters, which means the transformation process take considerable computing cost. • The transformation formulas for specific areas are based on estimation and inherently of errors, in Eq.(2.7); and therefore, positioning accuracy is degraded as such 11.

(18) Table 2.2: Parameters used in Eq.(2.1)–Eq.(2.15) : More specific parameters δX. 764.558m. δX. 361.229m. δZ. 178.374m. s. −0.2329 × 10−4. εX. 1.169855413 × 10−6. εY. −1.83986792 × 10−6. εZ. −9.822325179 × 10−7. K0. 0.9999. Ws. 250000m. λ0. 121◦. Sφ67. a67 · (A0 φ67 − A2 sin(2φ67 ) + A4 sin(4φ67 ) − A6 sin(6φ67 ) + A8 sin(8φ67 )). A0. 1 − 14 e267 −. A2. 3 2 (e 8 67. A4. 15 (e4 256 67. A6. 35 (e6 3072 67. A8. −135 8 e 131072 67. 3 4 e 64 67. + 14 e467 +. 5 6 e 256 67. 15 6 e 128 67. + 34 e667 − −. −. 41 8 e ) 32 67. −. 77 8 e ) 128 67. −. 175 8 e 16384 67. 455 8 e ) 4096 67. = 9.9832425786358 × 101. = 2.5146678795952 × 10−3. = 2.6391037362368 × 10−6. = 3.3889740139353 × 10−9. = −2.0687465476425 × 10−12. 12.

(19) GPS Package of Signal GPS Receiver. Transformation from Cartesian to Ellipsoidal in TWD97. φ84, λ 84, h 84. φ67, λ 67, h 67. Transformation from Ellipsoidal to Cartesian in WGS84. Transformation of 2o Transverse Mercatos Projection. X84, Y84, Z84 Datum Transformation from WGS84 to TWD67. E 6T7M 2. X67, Y67, Z67 Bursa-Wolf 7-parameter Transformation. N 6T7M 2. Figure 2.1: Level-wised transformation from WGS84 to TM2 in Taiwan errors are accumulated during the transformation process. • Traditional numeric methods [20, 21] can offer solutions for regression problems, but usually need to search the solution space exhaustively and take more computing time. • The more steps for transformations, the more inaccuracies. From the aspect of computational costs, the more steps for transformations may derive more power consumption; which is a serious problem when GPS is applied on a mobile device. The target of this study is to try to simplify the transformation equations in which the positioning accuracy still holds. A simpler form of the equations may be more efficient when used in instantaneous applications. Additionally, we expect the simplified transformation equations to be readable so that they can be easily encoded and integrated with GPS applications for further modifications. It is possible to derive solutions from these formulas by substituting the terms from one formula to another, but suffers from more computational costs and the accumulation of errors. Conceptually, the process in Figure 2.1 can be viewed as the composition of all transformation functions, which can be considered as a regressive function from source coordinates WGS84 to the target coordinate TWD67, as shown in Figure 2.2. If we can derive such formulas that perform. 13.

(20) GPS Package of Signal GPS Receiver. φ84, λ 84, h 84 Direct Regression Functions. E 6T7M 2 N 6T7M 2 Figure 2.2: Direct transformation from WGS84 to TM2 fewer levels of transformation, the inaccuracies and computational costs may be reduced. There are two types of techniques that can be applied for solving regression problems. Black-box methods mainly focus on the prediction accuracy. The mathematic insides of the regressive model are ignored since they are unnecessary or too complicated to be describable. On the other hands, white-box methods concerns both the prediction accuracy and the describability of the regressive models. In this study, we choose whitebox methods since the obtained regression model need to be describable and comparable with existing formulas.. 2.2. Related Work. Finding regressive models is essential in many applications, such as [22, 23, 24, 25]. Among these applications, black-box methods, such as artificial neural networks [26, 27] and support vector regression [28, 29], are widely employed. As mentioned before, blackbox methods are not dedicated to describing the relationships between the inputs and outputs. White-box methods offer the feasibility to describe the relationships between inputs and outputs in functional expressions. Traditional methods like linear regression [30, 31, 32] can serve as tools for finding regressive models. However, its accuracy and performance vary in different applications. Evolutionary-based methods like genetic algorithms and genetic programming are emerging methods for complex regressive problems, which usually gain satisfactory results under the control of sophisticated parameters. 14.

(21) In GPS applications, few studies devoted on using genetic-based methods for coordinate transformations. Below a brief review on related methods using genetic programming is given. [33] presented the implementation of symbolic regression based on genetic programming. In order to offer more possibilities for solving complicated problems in technological and scientific applications, Ferreira proposed “Gene Expression Programming (GEP)” which uses an array structure to represent the regressive tree [34]. [35] extended the mechanism of GEP to find classification rules in functional expression. [36] retained global information of the promising solutions during the genetic search process and used the information to control the operation of crossover. [37] improved the performance of genetic programming when searching for the symbolic regression via the analysis of dissimilarity and mating. [38] proposed a new selection strategy, called (µ + λ) - GP, of individuals for mutation, where λ descendants generated by µ parents will participate in the competition and be selected the best µ individuals as survivals.. 15.

(22) Chapter 3 The Design Concepts Symbolic regression can be implemented using the techniques of genetic programming. Genetic programming is a genetic-based algorithm which tries to find the combinations of operators and operands that can best describe the underlying problem. For a given problem, the input/output variables of the given problem are known as terminals. Meanwhile, a set of operators which possibly act as the components of the final expression √ is prescribed. For example, +, −, ×, ÷, sin, cos, n ·, ln, log, are typical operators in genetic programming. The tree-based implementation of genetic programming [39, 36] is employed, wherein the function to be found is represented as a binary tree with terminals as leaf nodes and operators and temporary variables as internal nodes. Each tree expresses a program or a function that describes the relationships between the input and out variables. A fitness function is usually defined as the least square error between the functional output and the expected output. The tree with lowest fitness value is retained. In this thesis, we follow the concept of tree-based implementation of genetic programming and extend the structure to meet the requirements of the underlying problem. First of all, we define the regression problem as follows.. 3.1. Problem Formulation. A regression problem is to find descriptive relationships between inputs and outputs so that the expected output corresponding to a new input is predictable. Formally, it can be described as follows. Suppose that D is a set of values appearing at the interval we are interested in and I = {(xi , yi ))|0 ≤ i ≤ n, xi ∈ D}. 16.

(23) is a set of paired data. The purpose of regression is to find a function g(x) so that yi ) ∼ = g(xi ) for all xi ∈ D. Usually, the error, εi , between f (xi ) and g(xi ) is defined as least square errors in the form of 1 εi = (yi − g(xi ))2 , ∀xi ∈ D. 2 To measure the effectiveness of g(x), an error function, E, is defined over all errors between f (x) and g(x), where E=. X. εi =. i. 1X (f (xi ) − g(xi ))2 , ∀xi ∈ D. 2 i. Regression is to find the best expressive function g(·) by minimizing E. One of the basic characteristics of applicable problems of regression is the existence of dependency between the paired data, i.e., (xi , f (xi )), for all i. Suppose G is the area we are interested in where GPS signals ~g in WGS84 are collected. Let IE = {(~ gi , TE (~ gi ))} and IN = {(~ gi , TN (~ gi ))}, ∀~ gi =hφ84 , λ84 , h84 i ∈ G, be the sets of paired data. TE (~ gi ) and TN (~ gi ) are the 2-D target eastings and northings collected from independent sources and correspond to g~i , respectively, in TM2. Let tE (~ gi ) and tN (~ gi ) be the values obtained by regression with respect to TE and TN , respectively. Similarly, let the easting-error εE (~ gi ) between TE (~ gi ) and tE (~ gi ) and the northing-error εN (~ gi ) between TN (~ gi ) and tN (~ gi ) be defined as square errors in the form of 1 (TE (~ gi ) − tE (~ gi ))2 , 2 1 εN (~ gi ) = (TN (~ gi ) − tN (~ gi ))2 , 2 εE (~ gi ) =. (3.1) (3.2). for all g~i ∈ G. Also, to measure the overall effectiveness of the regressive model on the easting points, the easting-error function EE is defined over all easting-error. The northing-error function EN is defined in a similar way. That is, for all g~i ∈ G, 1X (TE (~ gi ) − tE (~ gi ))2 , 2 g~i g~i X X 1 = εN (~ gi ) = (TN (~ gi ) − tN (~ gi ))2 . 2. EE = EN. X. εE (~ gi ) =. g~i. (3.3) (3.4). g~i. For easting-coordinates, minimizing EE may gain the best expressive function for tE (~ gi ). For northing-coordinates, a minimal EN is expected. Notably, EE and EN can be minimized independently if the results of eastings and northings are considered separately. If in an area the associated accuracy on both eastings and nothings are seriously asked, the area-error function EG can be defined as EG = wE × EE + wN × EN , 17. (3.5).

(24) where wE , wN ∈ [0, 1], wE + wN =1, are biases on EE and EN , respectively. Normally, wE =wN =0.5 for a region of similar lengths in both easting- and northing-sides. For areas where easting is narrow and northing is long, as in the area in this study, wE and wN should be tuned. Of course, the height above the ocean surface can be defined in the same way as Eqs.(3.1)–(3.5). However, it is not used in the underlying 2-D application. Next, a genetic-based method is proposed for deriving simpler formulas for direct transformation from WGS84 to TWD67 with respect to the minimization of EG .. 3.2. Tree-based Representation. As those in many other applications using genetic algorithms, the solution of the problem needs to be encoded as chromosome which is usually represented in bit-string or arrays, etc. However, for symbolic regression, unconstrained representation of chromosome may derive illegal or unreasonable solutions. For example, mutation on a string (2,+,5) which represents the addition of integer 2 and 5 may result in illegal expression (+,+,5). Since functional expressions can be parsed and analyzed in tree-structures, [39, 36] proposed binary tree representations for regression problems. Here, we define our extension. Let t be a binary tree, t=hV, Ei, where V is a set of vertex and E is a set of links. There are two disjointed sets of in V =VO ∪ VT . For convenience, we use O and T to denote the sets of operators and terminals. An element in Vo represents an internal node of t, addressing an operator used in the regression problem. An element in Vt represents a leaf which denotes a terminal used in the regression problem. In this study, we have two kinds of node or allele, terminals and operators. Terminals are φ84 , λ84 , h84 , and constant coefficients; operators are +, −, ×, ÷, sin, cos, etc. Additionally, each node is enumerated by its a location or locus Li in t by Eq.(3.6), noted by (Li , O or T ).  Lparent if node i is the left subtree of its parent   2 L Li = 2 parent + 1 if node i is the right subtree of its parent   1 if node i is root. (3.6). Figure 3.1 depicts some binary tree structures representing binary functional expressions. Note that, for simplicity, 1-ary operators, such as sin(x), exp(x), etc. are also represented as binary relations defined over two identical variables, e.g., sin(x, x) means sin(x) and exp(x, x) means exp(x), as shown in Figure 3.1(b)(c)(d). The definition of operators is shown in Table 3.1.. 18.

(25) ：+. ：. 1 2. 4. ：×. ：^. ：x. 1 +. ： sin. 3. 5. ：2. 3. 2. ：：. ：÷. ：x. 7 y 4 inv. 6. ：. ：. 1. ：. ：. 11. ：2. 14. ：2. 1. ：÷. ：y. 2y. ：2. 7. ：inv. ：y. 15. (d) legal, 2 × sin(exp(y)). ：×. 6. ：inv. ：+. ：×. 3 sqrt. 5. 14. 15. ：. ：x. ：y. ：exp. 6. 1. 2. 4. √. ：sin. ：inv. ：×. 7. (c) legal, ln(x − 2) −. 3. 2 2. 3. 6 inv. ：×. ：. ：sqrt. ：ln. 2. 10 x. 15. (b) legal, sin(x) + cos(xy). ：－. ：－. ：y. 14 x. 1. 5. ：×. 7. 13. (a) legal, x2 + ( y2 ) × y. ：inv. ：inv. ：2. 12 y. 4. 6. 5. ：cos. 3. 2. 7. ：×. 4. (e) illegal. ：sin. 5. ：y. 7. ：x. (f) illegal. Figure 3.1: Example of chromosome in tree-based representation, nodes are represented by locus:allele. inv means invalid node. 19.

(26) Table 3.1: Definition of operators in tree structure Binary. +. Function. Tree Expression. Expression. (Li , O or T ). x+y. (1, +) (2, x). −. x−y. ÷. ∧. sin. x×y. (1, ×). x×y. (2, x). x÷y. (1, ÷). x÷y. (2, x). xy. (1, ∧). xy. (2, x). Function. (Li , O or T ). sin(y). (1, sin) (2, invalid). cos. cos(y). (3, y) exp(y). (3, y) ln(y). (3, y). 20. ·. √. y. (3, y). (1, ln) (2, invalid). √. (3, y). (1, exp) (2, invalid). ln. (3, y). (1, cos) (2, invalid). exp. (3, y). Tree Expression. Expression. (3, y). (1, −) (2, x). ×. 1-ary. (3, y). √ (1, ·) (2, invalid). (3, y).

(27) 3.3. Defining Fitness. On the basis of the binary tree structures, this study designs a new fitness function and several improvements to improve the effectiveness of genetic-programming. Let t be a binary tree obtained in the process of symbolic regression for a region G. The fitness function f it(t) is defined over all errors between the target coordinates and the source coordinates, where f it(t) = w1 × f itv (t) + w2 × f its (t),. (3.7). ∀~ gi ∈ G. In Eq.(3.7), f itv (t) defines the functional errors based on Eq.(3.5), i.e., f itv (t) = EG = wE × EE + wN × EN .. (3.8). When f itv (t) is applied for tE and TE , wE =1 and wN =0; and conversely, wE =0 and wN =1 for tN and TN . In this study, wE =0.6 and wN =0.4 are used to reflect the geographic form of Taiwan. The indicator f itv (t) reflects if the expected regressive function is achieved. In Eq.(3.7), node(t) (3.9) 2d+1 − 1 defines the structural errors which considers the complexity of the binary t. In Eq.(3.9), f its (t) = 1 −. node(t) is the number of nodes in t and d is the depth of t. The indicator f its (t) controls the skewness of t. Usually, a complete and balanced binary tree is preferred and will have a smaller value of f its (t). Since considering on the correctness of the regression function and the structural complexity of the corresponding binary tree t is usually a trade-off, w1 and w2 are associated with f itv (t) and f its (t) in Eq.(3.7), respectively. By default, w1 , w2 ∈ [0, 1], w1 + w2 =1, are biases on f itv and f its , respectively. Moreover, we design one more fitness function for the purpose of computing power consumption. f ito (t) =. |O| X i=0. wi ×. Oi |I|. (3.10). where O is the set of operators using in the engine, |O| is the number of operators, and P Oi ∈ O, |O| i=0 Oi = |I|. The number of internal nodes in the tree t is denoted by |I|. wi , P wi =1, is defined by users who prefers final solutions with simpler operators, like + or √ −, or complicated operators, like · or ln. Consequently, the Eq.(3.7) will improved to the Eq.(3.11) f it(t) = w1 × f itv (t) + w2 × f its (t) + w3 × f ito (t),. (3.11). In the initial training phase, t is randomly generated and may be highly skewed; and, clearly, f it(t) is dominated by f itv (t), due to f its (t) ∈ [0, 1] but f itv (t) ∈ [0, ∞). After several iterations, f its (t) becomes effective as f itv (t) converges to [0, 1]. At the same time, f ito (t) will also effect f it(t) and press the engine to searching the solutions saving 21.

(28) more power. This is a reasonable design since it is more important to fast obtain an approximated solution rather than to make the binary trees complete.. 3.4. Genetic Operations. As in traditional GA applications, solutions are derived from parent chromosome by applying genetic operators, selection, crossover, mutation, and reproduction, etc. During the process of genetic programming, a number of trees are generated and retained for further processing. We define three pools for temporary results, Pev , Pco , and Pmt , where Pev keeps the surviving chromosome to be processed in the next generation, Pco retains the offspring after crossover, and Pmt is for the mutants. Let F be a set of trees on which the genetic operators will apply. In the following, we present the proposed operators.. 3.4.1. Crossover. Tree-based crossover is performed by swapping distinct sub-trees between the parents. In the pool of Pev , two main steps are performed in sequence. 1. Select from Pev two parents. 2. Determine which node in each parent tree are used for crossover. Suppose that ta and tb are two trees selected for crossover, where CPa and CPb are the nodes of ta and tb on which crossover is going to be applied. For convenience, CPa and CPb are referred to as ”crossover points” of ta and tb , respectively. Let t0 (t, p) represent the sub-tree of a tree t rooted by an internal node p of t. Tree-based crossover is to generate off-springs t0a and t0b from ta and tb , where • t0a =ta \t0 (ta , CPa ) ∪ t0 (tb , CPb ) and • t0b =tb \t0 (tb , CPb ) ∪ t0 (ta , CPa ). Figure 3.2 shows some examples of crossover. By defaults, the crossover point is randomly selected from a parent tree. We adapt tournament selection to select two randomly selected out of k trees; the one with better fitness value is retained for crossover. However, due to the effects of early-maturing, it is possible to derive diverse results and makes the process hard to converge. Several improvements are proposed for adaptively selecting parents and their corresponding crossover points according to their fitness. These will be presented in Section 3.4.3.. 22.

(29) ta +. ta′ +. tb + x. y. －. x. +. －. y. y. tb′. y. x. x. x. x. (a) CPa = 3, CPb = 3. ta. tb. ta′. tb′. +. +. +. y. x. y. －. x y. x. + x. －. x y. x. (b) CPa = 2, CPb = 1. y. ta′ +. tb +. ta +. x. －. x y. y. tb′ +. y. x. －. x x. (c) CPa = 3, CPb = 6. Figure 3.2: Example of crossover. 23. x.

(30) 3.4.2. Mutation. Tree-based mutation is performed by replacing either an internal node with another operator, a leaf node with a terminal other than the current one, or both in a subtree with new elements. As that in mutation, two main steps are performed in sequence. 1. Select from Pev a tree for mutation. 2. Determine which nodes in this tree will mutate. Suppose that ta is the tree selected for mutation and M Pa is a node of ta on which mutation is going to be applied. For convenience, such an M Pa is referred to as a ”mutation point” of ta . There may be more than one mutation points in a tree. Let rep(M Pa , e) denote that the operator/terminal defined in M Pa is replaced by a new operator/terminal e. Tree-based mutation is to modify a selected tree by one of the following operations. • V (ta )=V (ta )\{M Pa } ∪ {rep(M Pa , e)}, e ∈ O if M Pa ∈ O, • V (ta )=V (ta )\{M Pa } ∪ {rep(M Pa , e)}, e ∈ T if M Pa ∈ T , or • t0a =ta \t0 (ta , M Pa ) ∪ (t00 6= (ta , M Pa ) ). As the Figure 3.3 indicates, mutation provides another evolutionary behavior distinguishing from crossover. Mutation happens by considering whether the average fitness value of the all chromosomes found in the current population is improved, by comparing with the ones found in the past h populations. The probability of mutation defined as  1 − |fcit(j,j)−fcit(j−h,j)| if fc it(j, j) − fc it(j − h, j) > 0 0 fc it(j,j) R = random[0, 1.0] if fc it(j, j) − fc it(j − h, j) ≤ 0. (3.12). where fc it(i, j) is the average fitness value of all chromosomes found in the ith–jth population. If mutation happens in the jth population, a tree is randomly selected for mutation. The content of one or more nodes of the selected tree is replaced (operators by operators and terminals by terminals) or new sub-trees are generated and attached to it randomly. The same with crossover, selecting the mutation-point M P is random by default. Similarly, several improvements are proposed for adaptively selecting trees and their corresponding mutation points according to their fitness, which will be discussed in Section 3.4.3.. 24.

(31) ×. +. －. x y. －. x y. x. y. +. y. x. y. x. y. (b) Operators by operators,M P = 7. +. －. x. －. x. －. x. y. x. × －. x. －. x. (a) Terminals by terminals,M P = 1. +. +. +. y. x. y. (c) Two mutation points, M P = 1, M P = 7. ×. x ^ x. sin y. y. y. (d) Subtree Mutation, M P = 3. Figure 3.3: Examples of mutation. 3.4.3. Adaptive Selection. Generic Selection The operator “selection” is to select proper trees from F for a specific purpose, e.g., mutation or reproduction, and it can be applied at any place if needed. Let S(F, P r, n) denote the result of n trees by applying the selection operator on F with a probability method P r. For example, S(F, U, n) = {t1 , t2 , . . . , tn } is the result of n trees selected from F by applying uniform distribution U . In the initial stage of finding regressive models, generic selection operator works with uniform distribution for selecting trees or nodes. However, after some genetic iterations, such simple probability method may not derive significant result and need to be changed. In the following, we design several improvements on the generic selection operator to make it adaptive regarding to the environment.. Selecting Candidate Trees First of all, the basic of selective probability of a tree ti is defined as f it(ti ) P r(ti ) = 1 − P . i f it(ti ). 25. (3.13).

(32) P r(ti ) is directly considered when selecting trees for reproduction. The survivals will be saved in evolution pool, Pev . In the operation of reproduction, crossover, or mutation, trees are selected as candidates to be processed. After several number of iterations, the selecting method is changed to Eq.(3.13) by considering the goodness of each tree when selecting candidates. In order to further improve the effectiveness and efficiency of the proposed method, we present an adaptive mechanism for adjusting the crossover rate, Rco , and mutation rate, Rmt . The last k generations are observed and the tree with best fitness are retained in Pps on which we observe if the evolution process gets a decreasing performance. k is the given parameter and equal to the size of Pps . We devise an estimation function ftriple−D to evaluate such cases as follows. ftriple−D = w1 ∗ fdecreasing + w2 ∗ fdivergence + w3 ∗ fdemarcation where. 3 X. (3.14). wi = 1. The components of Eq.(3.14) are designed as below.. i=1. • Decreasing factor (fdecreasing ). This is for the case when the fitness decreases and rapidly converges in beginning but slows down or even becomes unchanged in the recent generations. fdecreasing = 1 −. M ax(f it(Pps )) − M in(f it(Pps )) M ax(f it(Pps )). (3.15). If fdecreasing is large, the regression engine becomes inefficient in searching solutions. In such cases, the probability of selecting trees for crossover or mutation should be increased. • Divergence factor (fdivergence ). This is the counter case of decreasing factor. That is, the fitness remains pretty high for a while but decreases and converges rapidly in the recent generations. Let gdivergence denote the number of generations wherein the best fitness remain unchanged. Intuitively, the larger gdivergence , the larger fdivergence is incurred. Therefore, fdivergence is defined as gdivegence , (3.16) fdivegence = k where k is the size of Pps . Also, if gdivergence is almost equal to k and fdivegence approaches to 1, the probability of selecting trees for crossover or mutation should be arisen. • Demarcation factor (fdemarcation ). This is to estimate if the average of the best fitness of the trees in the last k generation is stably converging but not early-maturing. Let Ldesired be a user-defined parameter. We define fdemarcation as  Ldesired 1 − Avg.f itness(Pps ) ≥ Ldesired ; Avg.f itness(Pps ) fdemarcation = 0 Avg.f itness(Pps ) < Ldesired 26. (3.17).

(33) After all, the adaptive selecting probability is revised by Eq.(3.18) with pre-defined thresholds adaptiveinc and adaptivedec as follows.   woriginal × R + wmove × (R + Rmove ) if ftriple−D ≥ adaptiveinc   0 R = woriginal × R + wmove × (R − Rmove ) if ftriple−D < adaptivedec    R otherwise. (3.18). The probability is adjusted according to the result of comparing Eq.(3.14) with adaptiveinc and adaptivedec . Biases woriginal and wmove are used to balance the effects of the old probability R or the new one adjusted. The new probability R0 is then applied for selecting tree for mutation or crossover.. Selecting Candidate Nodes For mutation or crossover, candidate trees are selected according to the process mentioned in Section 3.4.3. Suppose that a tree tc is selected as a candidate for crossover or mutation. Next step is to select from tc a node on which mutation or crossover will apply. Initially, all nodes are treated equally with uniform probability. After a number of generations, trees become complicated and possibly skewed. One of the requirements in the regression model is to simplicity and balanced. Therefore, we expect the trees obtained from mutation or crossover to be as simple and balanced as possible. We assume that a node vc in the skewed and complicated part of a subtree is likely to be useful in improving the fitness of tc if the genetic operators are applied on vc . The mechanism for selecting such a candidate node vc from tc adaptively is proposed as follows. • STEP 1. Separate the left- and right-subtrees of tc from the root node. We then L R L have two subtrees from tc , tR c and tc . The complexity of tc and tc can be calculated. by Eq.(3.9) individually. • STEP 2. In the subtree tcc with higher complexity, find the candidate node vc among the nodes of tcc . Each node in tcc is associated with a selecting probability fnode (vc , tcc ). fnode (vc , tcc ) =. 1 , L(vc , tcc ). (3.19). where L(vc , tcc ) denotes the level of vc in tcc . That is the candidate node is more likely to be close to the upper part of the candidate tree. This prevent candidate nodes from being located on the leaves and introducing higher complexity. • STEP 3. Repeat STEP 1&2 if more candidate nodes are necessary. Figure 3.4 show these steps. 27.

(34) tc. tc. or. tcR. tcL. tcL. tcR tcL > tcR. tcL < tcR. (a) Possible structure of skew subtrees.. tc tcL tcR. L=2. L=3 L=4. f node. L=5 L=6. CPc or MPc. L=7 R (b) Selecting CP or M P in ccc , when tL c > tc. Figure 3.4: Examples of selecting candidate node. 28.

(35) Table 3.2: All defined cases for modification Before. After. Before. After. Before. After. A+A. 2×A. A+0. A. A+1. unchanged. 0+A A−A. A×A. 0. A2. 1+A. A−0. A. A+0. 0−A. −A. 0+A. A×0. 0. A×1. 0×A A÷A. 3.4.4. 1. A. A. 1×A. A÷0. won’t happen. A÷1. A. 0÷A. 0. 1÷A. A−1. Tree Refinements. After several iterations, the trees with better fitness are retained in the pool Pev . The following methods are applied to prune trees to try to simplify their structures. • Constrains: Internal or leaf nodes are associated with some constrains for avoiding irrational representations. For example, ln(0), tan(π), zero-division, etc. Such irrational cases are rejected when they are produced in mutation or crossover. Also, terminals of 1-ary operators are restricted to be identical. • Modification on subtrees: In order to simplify the binary tree, a sub-tree is replaced by an equivalent, but simpler in structure, tree after several iterations. For instance, nodes representing x+x+x+x can be replaced by 4∗x and 3+5 by 8, etc. Figure 3.5 provides simple examples in graphic. However, only trees with better fitness values are considered to be modified. Furthermore, in Table 3.2, we show all modification cases in our engine. • Back tracking: If the best fitness value is retain stable in a long span, that means the engine is probably searching ineffectively and confronted by convergence of local optimal value. Back tracking will reload the survivals in the pass-n generation to create another opening for jump out the curve of local optimal. Clearly, we have a memory space, Pps , to store the survivals in the pass-n generations. The number of pass-n is bigger, there occurs more consumption of memory space. Users should think over the computing power of the machine to define the value of pass-n. 29.

(36) ＋＋. × A. ＋. 4. A. A. A. A (a) x + x + x + x replaced by 4 × x. × × + 3. × ×. A A. 8. A A. 5 (b) 3 + 5 replaced by 8.. × +. － A. × A. A. + 0. × A. A. A (c) A − A replaced by 0, and 0 + A replaced by A.. Figure 3.5: Examples of tree modification. 30. A. A.

(37) • Branch and bound: Branch and Bound is strict selection. The survival chromosome in Pev in generation i should be better than the best chromosome in Pbs such that the survival will pass the Branch and Bound and add to Pbs . The process provide the strict environment to press on better solution. But Branch and Bound won’t execute in every generation to avoid wasting too much computing resource and time cost as result of finding the best of the best. • Bound optimal local search: We believe that choosing the best-n chromosomes to combine with all operators will be born children better than parents. The matching process we called is Bound Optimal Local Search (BOLS). Every chromosome in Pbs will combine with primitive operators, therefore the computation order of BOLS is O(n × m) where n = max size of Pbs and m = number of operators. The process also will execute in a given circle, basically longer than 10 generations, rather than in every circle. • Constants: After modification, some chromosome will generate some nodes which’s value is the constant. If the leaf is constant and is calculated by addition or subtraction, then the engine will recombine the subtree. Moreover, the engine will sift the slide between maximum fitness value and minimum fitness value by addition to the negative average of fitness value of all chromosomes in the population. The process and the example are shown in the Figure 3.6.. 3.5. Termination Criteria. The termination criteria used are • Algorithm Converged: If there is no improvement in the search over the previous k generation, the convergence is exist. • Evolutionary time-out or over predetermined generations: The time of evolution and maximum generations are given by parameters of experiment. Of-course, more precise solutions need more computation time. On the contrary, too small generation number or time may not provide good solutions. • Expectant solutions found: The optimal solution may not be necessary, but it is good enough to satisfy the designer’s requirements.. 31.

(38) × × + 8. ×. − A A. ×. 3 A. A A. A. + 8 − 3. A A. Constant_of_modification = 5. (a) Modification first to get the constants node with addition and subtraction. ∑. Constant =. i∈Ninputs. . × (−1) + Constant_of_modification. Ninputs Min.  Constant  Fitness  Best chromosome   Information of tree  Root × × A. . [tE (gi ) −TE (gi )]. 0. Avg. Max. A A. A. (b) Adding the constants to the fitness value and shifting the fitness value. Figure 3.6: Constants in the engine. 32.

(39) 3.6. Processing Flow. Four main steps are applied iteratively in the training phase: • Generate an initial population of k binary trees, t1 , . . . , tk , by randomly composing the given operators and terminals as internal or leaf nodes. • Calculate the fitness value of ti , 1 ≤ i ≤ k, in the population. • Create a new population of k trees by reproducing the best existing trees or creating new trees by mutation or crossover (sexual reproduction of existing trees). • Satisfy the requirements and adopt the improvement strategies. • The best tree that appeared in any generation is designated as the result. With the definitions proposed above, the processing flow of our method is depicted in Figure 3.7.. 33.

(40) I E =<gi , TE (gi )>, I N =<gi , TN (gi )>, Operators Sampling Data. Calculation Fitness Value. tE (g i ), t N (gi ) εG ≤ Threshold. Simplifying Transformation Equations. Crossover Adaptive Probability Point selection with level-based selection. Mutation. Constraints Checking. Modification Back Tracking. Improvement?. Branch and Bound Bound Optimal Local Search. Selection of Best. Addition of Constant. Selection of Survivors. Reproduction. Figure 3.7: The proposed genetic-based process of finding regressive model. 34.

(41) Chapter 4 Implementation Issues This chapter describes the implementation of the ideas mentioned previously and discuss on the related problems when implementing such a system.. 4.1. Data Structures. We begin with defining the data structures representing tree-based chromosome. Treebased chromosome is directly implemented as a binary tree as shown in Figure 4.1. The components of such a tree-based structure are described in Table 4.1. A tree-based chromosome is a binary tree starting at Btree ptr. All chromosomes are using malloc to allocate the memory space and recycle after being used.. There are six main pools. being defined and used for various purposes, they are Pev , Pco , Pmt , Pbs , Pps , and Psv . All pools are defined as dynamic array, recording root-pointers, fitness value, etc., of Root. L_child 1 × 0 R_child. L_child 2 ＋ 0 R_child. L_child 4 NULL. x. 1 R_child NULL. L_child 5 NULL. L_child 3 － 0 R_child. y. 1 R_child. L_child 6. NULL. NULL. 2. 2 R_child NULL. L_child 7. y. 1 R_child. NULL. Figure 4.1: Tree-based chromosome: its data structure: (x + y) × (2 − y) 35. NULL.

(42) Table 4.1: The definitions of structures of a tree Name. Type. Comment. left. structure node. the pointer directs to the left subtree. right. structure node. the pointer directs to the right subtree. location. integer. the location of the node. value. double. the value of the node. type. integer. 0=operator, 1=argument, 2=constant, and -1=invalid. Pev. Evolutionary Pool 1 2 3 .. .. .. .. .. n. Tree chromosome. Pco. Crossover Pool. Pbs. Best Pool. Pmt. Mutation Pool. Pps. Pass-generation Pool. Psv. Survivors Pool.  Root  Information of fitness    Information of tree Constant   Other parameters. Left Child Location Value Type Right Child. Memory Space. Figure 4.2: Data structures of tree pools the trees under processing. All parameters are shown in Table 4.2 and Table 4.3. The data structures of the pools and their co-relations with functional trees are depicted in Figure 4.2.. 4.2. Implementing Genetic Mechanisms. There are four main procedures dealing with the temporary pools. We can get an overview of interaction among the procedures and tree pools in Figure 4.3. In this figure, long-dash lines are inputs, while the short-dash lines are the outputs. The following paragraphs will discuss on the details of the mechanisms and list their corresponding pseudo algorithms.. 36.

(43) Table 4.2: Parameters defined in tree pools Name root. Type Btree ptr. Comment the root of the tree. Information of fitness value fitness deviation. double. Eq.(3.8). fitness complexity. double. Eq.(3.9). fitness operators. double. Eq.(3.10). fitness. double. Eq.(3.11). fitness ori deviation. double. TTi (~ gi ) − tTi (~ gi ). Information of tree number node. integer. the number of nodes. number terminal. integer. the number of terminals. number addition. integer. the number of addition. number subtraction. integer. the number of subtraction. number multiplication. integer. the number of multiplication. number division. integer. the number of division. number others. integer. the number of others. level. integer. the number of level. generation. integer. the born generation of the chromosome. 37.

(44) Table 4.3: The definitions of structures of a tree Name constant. Type double. Comment the constants. Other parameters. 4.2.1. avg deviation. double. average deviation of population. max deviation. double. maximum deviation of population. min deviation. double. minimum deviation of population. kill. integer. 0 = survival, 1 = killed. execute crossover. integer. 0 = no crossover, 1 = crossover. execute mutation. integer. 0 = no mutation, 1 = mutation. The Initialization. The initialization procedure is to generate the primitive chromosomes which are the populations in the first generation. They are binary trees which randomly and legally combine input data pairs and operators. This can be described by the pseudocode presented in Appendix Algorithm 11. Note that, the Pev should be larger than |T |2 × |O| to have enough space to restore the primitives.. 4.2.2. The Selection Mechanism. The selection procedure is calculating fitness of trees retained in Pev , update all related information of trees, selecting the superiors to Pbs , and selecting the superiors to Pps .After the selection operation, only the trees with superior fitness survive, and continue to selection the survivals from Pbs and Pev in the procedure of selection of survivors. Improvement strategies like Branch and Bound, Selection-Best, Selection-Pass, and Selection-Survivorsare implemented as sub-modules in selection.The selection algorithms are described in pseudocode as shown in Appendix–Algorithm 12. The updating is sub-procedure of selection in Appendix Algorithm 13. The calculation of fitness needs the input data, denoted by IP, and is shown in Appendix–Algorithm 14. Branch and 38.

(45) Initialization. I. I. Selection. Fitness calculation Updating tree info. Selection Survivals O. O. Cross-over. Selection Parents Selection CPa CPb random Selection CPa CPb smart Cross-over 2p O. I. Mutation. Selection Mutants Selection MP random Selection MP smart Mutation node or subtree O. Pev Pco. I. Pmt Pbs Ppg. Improvement. Modification Constraint checking Back tracking Branch and bound Bound optimal local search Add constant O. Memory Space. Figure 4.3: Overview of symbolic regression engine. 39.

(46) Bound will be given a description later. In the selecting stage, the levels of trees and redundant copies of trees are examined. We devise two algorithms to deal with Level-Constraint in the pool Pev and the Check Redundancy in the pool Pbs . The max-level constraint help users to sieve out the undesired solutions. There is an especially noteworthy issue that the desired solution may not exist under the given level. The maximum level of a regressed function is constrained by users, which is checked by examining the information node of the tree structure. The algorithm of checking Level-constraint, Selection-Best, and Selection-Pass can be achieved by by simple if-then checking, so its pseudocode is ignored here. Evidently, checking redundant copies of trees in the Pbs pool takes a considerable amount of computational costs. We just examine the trees in the pool on their f it, f its , f itv , number node, number terminal, and level. If two trees with all information matched, redundancies of trees are found. The one in bad fitness is abandoned. The pseudocode of checking redundancy is presented in Appendix–Algorithm 15 and Appendix–Algorithm 16. Finally, we can select the survivors and put them into Psv for later operations. SelectionSurvivors is implements by following strategies, random selection, tournament selection, and truncation selection and uses them hybridly. There are two steps of the procedure. First step is selecting which pools, and second step is selecting the survivors from target pools. Target pools are Pbs , Pev , Pco , and Pmt . For simplification, we treat Pev , Pco , and Pmt as single larger pool. By given parameters, µ and λ, we define the proportion of Psv from the target pools respectively. In default setting, it is 50-to-50. If user sets the proportion manually, it’s so called truncation-selection. Then selection of survivors is implemented by random or tournament. We can define the number of players and choose the winner with best fitness value or roulette wheel strategy. The corresponding pseudocode is shown in Appendix–Algorithm 17, and the algorithm of tournament and one of roulette wheel will be presented later.. 4.2.3. The Crossover Mechanism. Crossover begins with selection of parents by rate Rco . If one of parent is selected for crossover, the engine will randomly find a spouse for it. Later in the generation, the selecting probability is adapted by the methods mentioned in Section 3.4.3. After completing the matching, we start to find the crossover points from each parent and start the crossover. The process is depicted by Figure 4.4 and Figure 4.5. Note that, we reserve twice spaces of Pev to be the size of Pco since there may be a lot of distinct trees generated with a higher Rco . The corresponding pseudocodes are presented in Algorithm 1 and Algorithm 2. In addition, the crossover is little different when CP is selected by level-based selection in Figure 4.6– 4.8 to prevent the trees to be skewer after crossover. 40.

(47) Random or Tournament. selecting a spouse. Pev 1 2 3 .. .. .. .. .. n. Copy to Pco(i+n). Pco 1 2 3 .. .. .. .. .. n .. .. .. .. .. .. .. .. 2n. Pev(i) Yes Rco. crossover? No Next chromosomes of Pev (a) Crossover. Pco 1 2 3 .. .. .. .. .. n .. .. .. .. .. .. .. .. 2n. (b) Crossover with two trees. Figure 4.4: Process steps of crossover The pseudocode is shown in Algorithm 3.. 4.2.4. The Mutation Mechanism. The mutation process begins with selection of a tree by probability method Rmt . If a tree is selected for mutation, the engine will randomly find a node to mutate. Later in the generation, the selecting probability is adapted by the methods mentioned in Section 3.4.3. The mutation operation is applied in two types: one-point mutation(terminal by terminal and operator by operator) and subtree mutation(terminal/operator by subtree). We select subtrees randomly or with tournament selection from Pev . Then we use the algorithm “level-based selection” for selection of M P to avoid choosing the small subtrees or leaves. In order to control the probability of replacing one-point or subtree 41.

(48) parent_ptr_1. parent_ptr_2. child_ptr_1. child_ptr_2 CPa. CPb. (a) Step.1. parent_ptr_1. parent_ptr_2. child_ptr_1. child_ptr_2 CPa. CPb. (b) Step.2. parent_ptr_1. parent_ptr_2. child_ptr_2. child_ptr_1 CPb. CPa. (c) Step.3. Figure 4.5: Steps of crossover-2p procedure. 42.

(49) Algorithm 1 Algorithm of crossover function crossover(Pev , Pco ) begin for each chromosome i of |Pev | do if random[0, 1.0] ≥ Rco then copy Pev (i) to Pco (i); choose a spouse to Pev (i) randomly or by tournament; copy spouse to Pco (i + |Pev |); end if ; end for; for each chromosome i of |Pco | do repeat selecting CPa of Pco (i) and CPb of Pco (i + |Pev |) randomly or with level-based selection; until CPa 6= 1 and CPb 6= 1; CALL crossover 2p(Pco (i), CPa , Pco (i + |Pev |), CPb ); end for; for each chromosome i of |Pco | do repairing the locus of trees in Pco ; end for; end;. 43.

(50) Algorithm 2 Algorithm of crossover with two trees 0. 0. sub-function crossover 2p(ta , CPa , tb , CPb ) begin create parent ptr 1, parent ptr 2, child ptr 1, and child ptr 2 ; create temp ptr ; 0. child ptr 1 points to the subtree of ta in locus=CPa ; 0. parent ptr 1 points to the subtree of ta in locus=b CP2 a c; 0. child ptr 2 points to the subtree of ta in locus=CPb ; 0. b parent ptr 2 points to the subtree of tb in locus=b CP c; 2. if child ptr 1 is left subtree then temp ptr point to left subtree of parent ptr 1 ; the left subtree of parent ptr 1 point to child ptr 2 ; end if ; if child ptr 1 is right subtree then temp ptr point to right subtree of parent ptr 1 ; the right subtree of parent ptr 1 point to child ptr 2 ; end if ; if child ptr 2 is left subtree then the left subtree of parent ptr 2 point to temp ptr ; end if ; if child ptr 2 is right subtree then the right subtree of parent ptr 2 point to temp ptr ; end if recycle(parent ptr 1, parent ptr 2, child ptr 1, child ptr 2, and temp ptr ); end;. 44.

(51) or. Level (ta ) ≥ Level given. taL. taR. CPa. taR. taL Level-Based Selection. Level (tb ) < Levelgiven. tbL. CPb. taR. Random Selection. Figure 4.6: CPa is selected by level-based selection, but CPb isn’t. Level (ta ) < Levelgiven. taL. CPa. taR. Random Selection. or. Level (tb ) ≥ Levelgiven. tbR. tbL. CPb. tbR. tbL Level-Based Selection. Figure 4.7: CPb is selected by level-based selection, but CPa isn’t. 45.

(52) Algorithm 3 Algorithm of selection crossover-point with level-based selection sub-function crossover level based selection(ta , CPa , tb , CPb ) begin if Level(ta ) ≥ Levelgiven and Level(tb ) < Levelgiven then CPa ← level-based point selection(ta ); CPb ← random selection; end if ; if Level(ta ) < Levelgiven and Level(tb ) ≥ Levelgiven then CPa ← random selection; CPb ← level-based point selection(tb ); end if ; if Level(ta ) ≥ Levelgiven and Level(tb ) ≥ Levelgiven then CPa ← level-based point selection(ta ); if CPa is the node of left subtree of ta then CPb ←level-based point selection(right subtree(tb )); end if ; if CPa is the node of right subtree of ta then CPb ← level-based point selection(left subtree(tb )); end if ; end if ; end;. 46.

(53) or. Level (ta ) ≥ Levelgiven. CPa. t. taL. taR. L a. taR. Level-Based Selection. or. Level (tb ) ≥ Levelgiven. CPb. tbL. tbR. tbL. tbR. Level-Based Selection. Figure 4.8: CPa and CPb are selected by level-based selection in locus M P , a parameter Rt p is re-defined. This process is depicted by Figure 4.9 and Figure 4.10 and described by the pseudocodes in Algorithm 4. Its sub-procedures are presented in Appendix–Algorithm 18 and Appendix–Algorithm 19.. 4.2.5. The Improvement Mechanism. • Back Tracking Mechanism: If the trees retain in Pbs for a long time, it seems that the engine falls into local optima. We can use the Eq.(3.14) to evaluate the status of convergency. If it happens, we pick up from Pbs some tree by probability calculated by Eq.(3.13) and put them to Pev . Backing Tracking provides a possible way to ruffle the current evolutionary pool and may give a chance to jump out from the local optima. The pseudocode is presented in Appendix–Algorithm 20. • Branch and Bound Mechanism: Branch and Bound can look as a stricter selection, which asks that the chromosomes in Pev in every generation must be better than the best chromosomes in Pbs . If Branch and Bound doesn’t execute, the selection best strategy will work instead. The selection best strategy lets Pbs restore the best chromosomes from Pev which are better than the worst ones in Pbs in every generation. The pseudocode is presented 47.

(54) Pev 1 2 3 .. .. .. .. .. n. Pmt 1 2 3 .. .. .. .. .. n. Pev(i) Yes Rmt. mutation? No Next chromosomes of Pev (a) Selecting chromosomes to be mutated. Pev 1 2 3 .. .. .. .. .. n Pmt. 1 2 3 .. .. .. .. .. n. Random or tournament candidate i subtree Pev(i). mutation type?. selection cutting point randomly or level-based selection subtree mutation Rt_p. point mutation. 1 2 3 .. .. .. .. .. n. 1 2 3 .. .. .. .. .. n. Setargument. Setoperator (b) Mutation with two types. Figure 4.9: Process steps of mutation. 48.

(55) parent_ptr_1. parent_ptr_1. child_ptr_1. child_ptr_1. ＋－ × ÷ … ln √. MP. random. (a) Operators by operators. parent_ptr_1. parent_ptr_1. child_ptr_1. child_ptr_1. MP. x y … z. constant. MP. random. (b) Terminal by terminal. parent_ptr_1. parent_ptr_1. Pev t1 t2 t3 … … tn. child_ptr_1 MP. random or tournament. (c) Replace by subtrees. Figure 4.10: Examples of mutation. 49.

(56) Algorithm 4 Algorithm of mutation function mutation(Pev , Pmt , T , O) begin for each chromosome i of |Pev | do if random[0, 1.0] ≥ Rmt then copy Pev (i) to Pmt (i); end if ; end for; for each chromosome i of |Pmt | do selecting M P of Pmt (i) randomly or with level-based selection; if random[0, 1.0] ≥ Rt p then CALL mutation subtree(Pmt (i), M P , Pev ); end if ; else CALL mutation point(Pmt (i), M P , T , O); end else; end for; for each chromosome i of |Pmt | do repairing the locus of trees in Pmt ; end for; end;. 50.