170 likes | 365 Views
Genetic Programming: Finding Approximate Analytic Solution to First-Order Ordinary Differential Equations. Differential Equations. Initial Value Problem: Find y=f(x) if y’ = f( x,y ) at y(a)= b (1) Closed-form solution: explicit formula -y’ = y at y(0)=1 (Separable)
E N D
Genetic Programming:Finding Approximate Analytic Solution to First-Order Ordinary Differential Equations
Differential Equations Initial Value Problem: Find y=f(x) if y’ = f(x,y) at y(a)=b (1)Closed-form solution: explicit formula -y’ = y at y(0)=1 (Separable) Ans: y = e^x (2) Numerical solution: y’ = (x*y)/(x+y) + y^2 at y(0)=1 Can’t solve explicitly… So we’ll use numerical methods: (a) Dormand-Prince (b) ode45: returns [x,y] values
Regression / Curve-fitting • Given data points, find best fitting curve • Normal regression finds coefficients of some predefined function (ie linear regression) • Symbolic Regression finds function and coefficients • Since normal requires human intuition, we’ll use symbolic regression using genetic programming
Genetic Programming Basics Tree Structure Function Nodes from a Function Set Terminal Nodes (numbers, variables) from a Terminal Set Create Random Function Trees Making the function: y = x*sin(3.63/(2+x))
Function Tree = Individual • Use Polish Notation • Operator comes first. Ex: 5 + x – (x * 8)) + 5 – x * x 8 Quiz: + 3 * 6 – 2 1 = ? • Syntax-preserving: binary or unary functions? • Binary: takes two nodes • Unary: takes one node
Fitness Evaluation Evaluate error at each point Fitness= -Total Sum of Error If we want greater accuracy, we will want to evaluate more points in domain We may also want to limit the maximum error between any sample point and the function
Reproduction • Random tournament to select who mates • Randomly select two individuals in population • The one with the best fit is parent 1 • Randomly select crossover or mutation • If mutate, mutate parent 1 • Else find parent 2 and crossover • Repeat reproduction until Population size is filled
Reproductive Operator: Crossover • Crossover (GP) Crossover (DNA analogy) Parent 1: Parent 2: (a+b*sqrt(b)) / sqrt(b+b) a*(((b/b)/a) * sqrt(b+b)) Child: a+(b*(b/b)/a) / sqrt(b+b) ( Throw away second child)
Reproductive Operator: Mutation Randomly select node in tree to mutate Randomly select a node of the same type Mutate node Mutation gets us out of pitfalls
Parameters • Determine Number of Data Points in Domain (Matlab) • Step size: 0.001 • Determine Function Set • F = {sin, cos, e^, /, *, +, -} • Determine Terminal Set • T = {x, pi/2, 100 random numbers between -5 and 5} • Create a Random Population of Function Trees • 10,000 indiv. w/ initial depth of 5 • Set Minimum Tolerance (error) • Tolerance: 0.01 • Set Maximum Num of Generations • Max num of generations: 600 • Set Crossover and Mutation rate • Crossover: 80% Mutation: 1-Crossover = 20%
Results: Finding a Closed-form Solution • y’=-y*cos(x) Solvable Ans: y=1 / e^(sin(x)) • Using Genetic Programming • Pop size: 1000 • Domain [0 , 5] • Problem • Protected division
Results: Finding an Approximation • y’=cos(2*x*y^2) • No closed form solution • Polyfit fails in matlab
Results: y’ = cos(2*x*y^2) …the explicit solution is…
Program Evaluation …the explicit solution is…
The Approximate Solution y = e^( ((((x / cos( sin( ((cos( -1.3977331148909333) / (x * ((x - (sin( e^( 0.5839073910575703)) / ((0.5839073910575703 + (((cos( -1.0768808474282068) / (((x - e^( (sin( (x + (-2.1183334622298835 + cos( -1.3977331148909333)))) * (sin( (x + 4.404029716618364)) / (x + x ))))) / -3.9294994016950726) * cos( cos( cos( ((x + -2.1183334622298835) / (((-3.488290785681106 / 4.429086745576107) - (4.64613256685335 + (2.032700600098037 / (cos( x ) / (3.951688988841374 / cos( sin( (((-1.6691938170386078 / cos( ((2.032700600098037 / -1.3977331148909333) / cos( cos( sin( (0.5839073910575703 + (0.5839073910575703 + (cos( 0.21588984462299265) / cos( 0.5839073910575703)))))))))) + 0.5839073910575703) + (0.5839073910575703 + (cos( 0.21588984462299265) / cos( 0.5839073910575703))))))))))) + -2.1183334622298835))))))) / -1.3977331148909333) + (cos( ((0.5839073910575703 / (0.5839073910575703 * 0.5839073910575703)) * (sin( (x + e^( 0.5839073910575703))) / (x - cos( e^( 0.5839073910575703)))))) / cos( 0.5839073910575703)))) / (3.951688988841374 / cos( sin( (0.5839073910575703 + (0.5839073910575703 + (cos( 0.21588984462299265) / cos( 0.5839073910575703)))))))))) * -1.3977331148909333))) / (sin( -4.09932580082895) / -1.3977331148909333))))) / cos( sin( ((cos( -1.3977331148909333) / (x * ((x - (sin( e^( 0.5839073910575703)) / ((0.5839073910575703 + (((cos( -1.0768808474282068) / (((0.5839073910575703 - e^( (sin( (x + -2.1183334622298835)) * (sin( (x + 4.404029716618364)) / (x + x ))))) / -3.9294994016950726) * cos( cos( cos( 0.5839073910575703))))) / -1.3977331148909333) + (cos( (e^( 0.5839073910575703) * (sin( (x + e^( 0.5839073910575703))) / (x - cos( e^( 0.5839073910575703)))))) / cos( 0.5839073910575703)))) / (3.951688988841374 / cos( sin( (0.5839073910575703 + (0.5839073910575703 + (x / cos( 0.5839073910575703)))))))))) * -1.3977331148909333))) / (sin( -4.09932580082895) / -1.3977331148909333))))) - e^( (sin( (cos( (cos( (((0.5839073910575703 / (0.5839073910575703 * 0.5839073910575703)) * (sin( (x + sin( (((2.032700600098037 / -1.3977331148909333) / 0.5839073910575703) + -2.1183334622298835)))) / (x - cos( e^( 0.5839073910575703))))) + -1.6691938170386078)) / (((e^( x ) + sin( e^( (x / cos( cos( (x / e^( 0.5839073910575703)))))))) - ((x / cos( sin( (((0.5839073910575703 - cos( -1.3977331148909333)) / (x * (3.141592653589793 * -1.3977331148909333))) / (x / -1.3977331148909333))))) / -1.3977331148909333)) + ((0.5839073910575703 / ((x / cos( cos( x ))) - (0.5839073910575703 * cos( cos( (0.5839073910575703 / (((-3.9294994016950726 * sin( x )) - cos( (0.5839073910575703 - (0.5839073910575703 * cos( (x / 0.5839073910575703)))))) + -2.1183334622298835))))))) - cos( -1.3977331148909333))))) + (x + sin( (((2.032700600098037 / -1.3977331148909333) / sin( 0.5839073910575703)) + -2.1183334622298835))))) * sin( (cos( (x / (cos( (x - (((cos( -1.0768808474282068) / (((x - e^( (sin( (x + (x - ((x - (sin( e^( 3.9169117826293363)) / ((0.5839073910575703 + (((cos( -1.0768808474282068) / (((3.141592653589793 - e^( (sin( (x - (sin( (x - 0.5839073910575703)) + 0.5839073910575703))) * (sin( (x + 4.404029716618364)) / (x / 4.5366847258566025))))) - -3.9294994016950726) + cos( cos( cos( ((x - -2.1183334622298835) / (((-3.488290785681106 / -3.423767380013757) / (4.496363738665307 + (2.032700600098037 + (cos( -2.9602495684565455) / (3.951688988841374 / (x - cos( e^( 0.5839073910575703)))))))) + -2.1183334622298835))))))) * -1.3977331148909333) * (cos( ((0.5839073910575703 / (0.5839073910575703 * 0.5839073910575703)) * (sin( (-1.7908824571411053 + e^( 0.5839073910575703))) / (x - cos( e^( x )))))) + e^( 0.5839073910575703)))) / (3.951688988841374 / cos( sin( (0.5839073910575703 + (-2.8354922153845505 + (cos( 0.21588984462299265) / cos( 0.5839073910575703)))))))))) * -1.3977331148909333)))) * sin( (cos( cos( 0.5839073910575703)) + (sin( 0.5839073910575703) / (x - ((((2.032700600098037 / -1.3977331148909333) / sin( ((sin( (cos( (x - (x + e^( 0.5839073910575703)))) + sin( (x + cos( sin( e^( (0.5839073910575703 * (sin( (x + 4.404029716618364)) / (x + x )))))))))) / cos( x )) / x ))) + -2.1183334622298835) + 0.5839073910575703)))))))) / -3.9294994016950726) * cos( 0.5839073910575703))) / -1.3977331148909333) / ((x - (cos( 0.5839073910575703) + (sin( e^( 0.5839073910575703)) / (((x - (sin( e^( (((0.5839073910575703 / ((x / cos( cos( x ))) - (0.5839073910575703 * cos( (cos( 4.404029716618364) / (x + cos( (-1.6691938170386078 / x )))))))) - cos( (sin( (x - cos( ((0.5839073910575703 / (0.5839073910575703 * (sin( 0.5839073910575703) * cos( 0.5839073910575703)))) * (0.5839073910575703 / (0.5839073910575703 * 0.5839073910575703)))))) * -1.3977331148909333))) / -4.09932580082895))) / (4.5366847258566025 / (3.951688988841374 / cos( -1.3977331148909333))))) * -1.3977331148909333) - cos( e^( 0.5839073910575703)))))) / cos( cos( (x - cos( e^( 0.5839073910575703))))))))) + 3.951688988841374))) + (sin( e^( 0.5839073910575703)) / (x - cos( e^( 0.5839073910575703))))))))) / -3.9294994016950726)) Notsimplified (MaximanorMatlabcouldsimplify) Fitnesspenaltyforlength
Sources Adapted Tiny_GP Java code RiccardoPoli. Found at http://www.gp-field-guide.org.uk/ Burgess, Glenn. “Finding Approximate Analytic Solutions to Differential Equations” commissioned by Department of Defense (Australia) in 1999 Koza, John. Genetic Programming RiccardoPoli