Computer Algebra with Mathematica:  Final Project

Genetic Algorithm for Optimization

Kevin C. Connolly
kcconnolly@aol.com
Instructor:  K. Levasseur

This project was completed during the Summer 1999 session.  I've made a few cosmetic changes to the notebook that Kevin summitted. - K. Levasseur

Note:  The first initialization cell appears in section II.C1 below.  Once all initialization cells have been evaluated, you can go to section IV.E, "Putting it all together," to see a sample run of the algorithm, and to launch further runs of your own.

I:  Background

The problem to be investigated in this project arises from a study of calling card usage done for a major telecommunications carrier.  Respondents were asked several questions about their current calling card usage (including mean monthly number of calling card calls, which in this sample turned out to be 7, and minutes, which turned out to be 24), and then were asked to rank in order of preference 16 different calling card offers, each of which had some combination of values from five key variables:  brand (Carrier 1 vs. Carrier 2); access method (1-800 access vs. 0+ access); monthly fee (0, $2.00/month to cover first 9 calling card minutes, $4.00/month for 20 minutes, $8.00/month for 45 minutes); per-call surcharge (0, $.25, $.50, $1.00); and per-minute rate (0, .10, .25, .45).  The configurations of the different offers were chosen to facilitate conjoint analysis.  The analysis then yielded utilities for each level of each variable, enabling us to see how much weight each variable carries in customers' decision-making process, and to see the relative impact of the different levels of each variable.

Each respondent was also asked to indicate how many of the 16 calling card offers they would consider acceptable for usage in their household.  Based on all the information collected, we were able to determine an average revenue, and then an effective revenue for each plan.

Example:  A plan with no monthly fee, a per-call surcharge of $.50, and a per-minute rate of $.25 would yield an average revenue of $9.50.  (7 calls @ $.50/call: $3.50  + 24 minutes @ $.25 per minute : $6.00 = $9.50).  If 40% of respondents indicated that this offer would be acceptable for use in their household, then the effective revenue would be $3.80 (average revenue of $9.50 * .40  = effective revenue of $3.80).

The goal of this project is to use a genetic algorithm to "evolve" a plan that will maximize effective revenue.  The project will have two major phases.  First, it will be necessary to develop a regression model from the collected data so that we can estimate customer preferences with regard to combinations of values that were not explicitly tested (and, using interpolation, with regard to values that were not tested at all).  This will enable us to compute an effective revenue figure for any offer, and thus to establish its "evolutionary fitness," with the fittest offers being those which yield the highest effective revenues.

The second phase will be the creation and execution of the algorithm itself.  We will begin with a randomly generated population of offers.  At each generation, we will test each member of the population for fitness (i.e., effective revenue), select according to fitness for reproduction, and apply the evolutionary mechanisms of crossover and mutation.  We will allow "evolution" to take place for fifty generations, and chart the effective revenue yielded by the fittest member of each generation.

II:  Regression Model

The utilities from the conjoint model are normed to 0, so that positive numbers indicate options which respondents tend to find attractive, while negative numbers are attached to those values which incline respondents away from particular offers.  The higher the number, the more appealing it is to respondents.

A & B.  Brand & Access Method

While these two variables do not affect average revenue, they do affect customer preference, and thus effective revenue.  Unlike the others, they are treated as discrete rather than continuous variables.  Carrier 1 has a utility of 0.3; Carrier 2 has a utility of -0.3.  Dialing 0 + number for access has a utility of 0.1; requiring the customer to dial a 1-800 number for access has a utility of -0.1

C.  Monthly Fee

This is the most complicated variable, because it requires two regressions:  one to fit the amount of the fee to the number of minutes it covers, and one to fit the amount of the fee to the utilities.  We'll take them in that order.

C1.  Fee for Minutes regression

<< Graphics`Colors`

minutesforfee  =  {{0, 0},  {2, 9},  {4,  20},  {8,  
            45}};

dataplot1 = 
      ListPlot[minutesforfee,  PlotStyle -> {Blue,  PointSize[ .02]}]

minutesforfeeEq  = Fit[minutesforfee,  {1, x, x^2},  x]

dataplot1a  = 
      Plot[minutesforfeeEq,  {x,  0,  8},  PlotStyle -> Red]

Show[dataplot1,  dataplot1a]

C2.  Utilities for Fee Regression

As can be seen below, customers strongly preferred not to have to pay a regular monthly fee for calling card usage, assigning the lack of a monthly fee a utility of 2.6.

monthlyfeeutilities = {{0, 
            2.6},  {2,   .2},  {4,  - .4},  {8,  -2.4}};

monthlyfeeutilitiesEq  = 
      Fit[monthlyfeeutilities,  {1, x, x^2},  x]

dataplot2  = 
      ListPlot[monthlyfeeutilities,  PlotStyle -> {Blue,  PointSize[ .02]}]

dataplot2a  = 
      Plot[monthlyfeeutilitiesEq,  {x, 0, 8},  PlotStyle -> Red]

Show[dataplot2, dataplot2a]

D.  Per-call surcharge

Customers prefer the surcharge to be as low as possible.  However, the fact that they assign the lack of a surcharge a lower utility than the lack of a monthly fee indicates that they feel more strongly about the monthly fee, and are willing to make the trade-off in order to avoid a regular commitment.

surchargeutility = {{0,  
            1},  { .25,   .7},  { .5,  - .4},  {1,  -1.3}};

surchargeutilityEq  = 
      Fit[surchargeutility,  {1,  x,  x^2},  x]

dataplot3  =  
      ListPlot[surchargeutility,  PlotStyle -> {Blue,  PointSize[ .02]}]

dataplot3a  = 
      Plot[surchargeutilityEq,  {x,  0,  1},  PlotStyle -> Red]

Show[dataplot3, dataplot3a]

A note on quadratic regressions

I chose to use quadratic regressions as a happy medium between what I perceived to be disadvantages on either side.  On the one hand, quadratic regressions provide a closer fit than do linear regressions.  On the other hand, while cubic or higher-order regressions do  provide exact fits at the given data points, they also have a tendency in my experience to overfit the data in such a way as to create a misleading shape to the curve.  Fitting the surcharge utility data to a cubic regression provides one example:

surchargeutilitycubicEq  =  
      Fit[surchargeutility,  {1,  x,  x^2,  x^3},  x]

dataplot3cubic  =  
      Plot[surchargeutilitycubicEq,  {x, 0, 1}]

The "false" dip at about x=.85 is what I'm trying to avoid--I don't think it fairly represents customer preference.  Except for your lecture notes on approx.nb, however, I'm largely self-taught in this regards, so if there are any considerations I'm neglecting, just let me know.

I guess that the main thing is that your approximating function i s decreasing (or maybe increasing in some cases)  I can't say that I'm familiar enough with regression theory to say whether a certain degree is better than another.

E.  Per-minute rate

Again, customers favor the lowest possible rates.  The fact that the range of utilities is broader here than for the per-call surcharge indicates that customers are willing to accept higher surcharges to get lower per-minute rates.

perminuteutility = {{0,  
            2.5},  { .1,   .8},  { .25,  -1},  { .45,  -2.3}};

perminuteutilityEq  = 
      Fit[perminuteutility,  {1, x, x^2},  x]

dataplot4  = 
      ListPlot[perminuteutility,  PlotStyle -> {Blue,  PointSize[ .02]}]

dataplot4a  = 
      Plot[perminuteutilityEq,  {x,  0,   .5},  PlotStyle -> Red]

Show[dataplot4, dataplot4a]

F. Connecting utilities with customer preference

As mentioned in the background section above, 16 offers were explicitly tested with respondents, and respondents were asked to indicate how many of the offers they considered acceptable for use in their household.  By pairing the total utility (brand + access + fee + surcharge + rate) of each offer with the percentage of customers who found that offer acceptable, we're able to develop a regression that will allow us to estimate the acceptability of those offers we have not explicitly tested.

offerutilities = {{4.2,   .72},  {2.5,   .57},  {1.3,  .41},  {1.0,   .37},  {1.0,   .37},  {- .7,  .32},  {- .2,   .27},  {1.2,   .34},  {2.9,   .45},  {- .9,   .22},  {1.9,   .45},  {1.3,   .36},  {-2.7,   .18},  {-4.2,   .11},  {-1.6,   .17},  {-3.1,   .09}};

offerutilitiesEq  = 
      Fit[offerutilities,  {1,  x,  x^2},  x]

dataplot5  = 
      ListPlot[offerutilities,  PlotStyle -> {Blue,  PointSize[ .02]}]

dataplot5a  =  
      Plot[offerutilitiesEq,  {x,  -4.2,  4.2},  PlotStyle -> Red]

Show[dataplot5, dataplot5a]

III:  Measuring Fitness

The variable we want to maximize from the offers is effective revenue, defined as the revenue which each offer would produce from the average consumer, multiplied by the percentage of consumers who find that offer acceptable.  We thus need to create functions which will allow us to compute this for each member of the "population," defined as an offer.

A.  Creating offers as members of the population

Each offer, as we've seen, has to have five components.  We'll take a shortcut in the definition of the brand and access method, allowing them to be represented immediately by their utilities.  We'll let the monthly fee be defined as a real number between 0 and 8; the per-call surcharge as a real number between 0 and 1, and the per-minute rate as a real number between 0 and .45   Thus

member  :=  { .3 * (-1)^
            Random[Integer,  {1, 2}],   .1 *  (-1)^
            Random[Integer,  {1, 2}],  Random[Real,  {0, 8}],  
        Random[Real],  Random[Real,  {0,  .45}]}

Creating a population of five lets us see that we get the output we expect

testpopulation =  Table[member,  {5}]

B.  Computing the revenue from an offer

As we've seen, the average calling card customer uses 24 minutes per month.  If, then, the customer's monthly fee pays for 24 or more minutes per month, then the monthly fee will be the sole source of revenue from that customer.  We can use the "minutes for fee" regression developed above to determine what level of monthly fee will buy the customer 24 minutes:

revcutoff  = 
      Max[x  /. Solve[- .0818182  +  
                4.32045 x  +   .164773 x^2  ==  24]]

We can discard the negative solution:  if the customer's monthly fee is greater than $4.72, then the average revenue from that customer will simply be equal to the monthly fee.  Since the monthly fee appears as the third element of each member of the population, we have:

avgrevenue[y_List  /;  y[[3]]  > revcutoff]  :=  
      y[[3]]

If the monthly fee is less than or equal to $4.72, the problem becomes slightly more complicated.  We have to determine how many calling card calls and how many calling card minutes are not covered by the monthly fee, using the estimate from the original data set that each call lasts approximately 24/7 minutes.  Then the revenue from the per-call surcharge and the per-minute rate for the portion not covered by the monthly fee must be computed.  This gives us the following:

avgrevenue[y_List ] :=  
      Part[y, 3]  +  (Part[y, 
              4] * (24  -  (minutesforfeeEq  /.  
                    x -> Part[y, 3])) * 
            (7/24))  +  (Part[y, 
              5] *(24  -  (minutesforfeeEq  /.  
                    x -> Part[y, 3])))

With this equation, we would expect an offer with no monthly fee, a $0.50 per-call surcharge, and a $0.10 per-minute rate to give us revenue of $5.90:  (7 calls at $.50 per call + 24 minutes at $.10 per minute = $5.90).

testcall = {- .3,  - .1,  0,   .5,   .1}

avgrevenue[testcall]

As can be seen, we get something very close to this, the slight discrepancy being due to the small constant term (-.0818182) in the minutes-for-fee regression.  Now that we can compute the average revenue from any offer, we are thus halfway to being able to compute effective revenue--and thus evolutionary fitness.  The second half involves the computation of the acceptability of each offer.

C.  Computing the acceptability of an offer

The various regression equations for utilities we developed in section II above now come into play.  First, we generate the total utility of an offer.  Here it's important to keep in mind that the brand and access method of each offer are already expressed as utilities in the representation of each offer.  Thus:

totalutility[y_List]  :=  
      Part[y, 1]  +  
        Part[y, 2]  +  (monthlyfeeutilitiesEq  /.  
            x ->  Part[y, 3])  +  (surchargeutilityEq  /.  
            x -> Part[y, 4])  +  (perminuteutilityEq  /.  
            x -> Part[y, 5])

We've also already developed the equation (in section II.F above) that will allow us to determine what percentage of the population will find an offer acceptable, given its total utility.  Thus:

acceptability[y_List]  :=  
      offerutilitiesEq  /.  x -> totalutility[y]

D.  Computing fitness

The fitness of an offer--its effective revenue--can be determined simply by the product of its average revenue and its acceptability.

effectiverevenue[x_List]  := 
      avgrevenue[x] acceptability[x]

testpopulation  =  Table[member,  {5}]

Map[effectiverevenue,  testpopulation]

Mapping effectiverevenue over a small test population gives us confidence that we are getting correct results.  The offer with the highest effective revenue is the one with the lowest monthly fee (and we've seen that avoiding the payment of a monthly fee is the most important consideration to customers).  The ones with the lowest effective revenue are those either with very high monthly fees (as in the last member of the population), or moderately high monthly fees combined with higher per-calls of other charges.

IV.  The Genetic Algorithm

A.  Introduction

As initially developed by John Holland in Adaptation in Natural and Artificial Systems (2nd ed--Cambridge, MA:  MIT Press, 1992), and as expounded by Melanie Mitchell in An Introduction to Genetic Algorithms (Cambridge, MA:  MIT Press, 1996), there are three major steps in the implementation of a genetic algorithm.  First is the decision regarding encoding the various members of the population under consideration so that "evolution" can work on them.  Second is the decision concerning selection--which members of the population in a given generation will have the opportunity to pass their "genes" on to the next generation.  And third is the decision of how genetic operators (such as crossover and mutation) will be allowed to work on members of a population within a given generation to diversify and strengthen the gene pool.

B.  Encoding

We've already tackled the problem of encoding in section IIIA above.  To review, each offer will be expressed as a list with five elements.  Brand and access method will be expressed in terms of their utilities; monthly fee, per-call surcharge, and per-minute rate will be expressed as decimal representations of their amounts.  Thus:

■Offer already defined, encoded, and initialized as member above

samplemember  =  Table[member,  {1}]

The sample member of the population encodes an offer from Carrier 2, offering 0+ access to an operator, with a monthly fee of approximately $1.96.  Calls not covered within the monthly fee will be billed with a surcharge just over 48 cents per call, and a per-minute rate of almost 13 cents.

C.  Selection

In An Introduction to Genetic Algorithms (henceforth IGA), Mitchell outlines several methods of selecting which members of a population can pass into the next population, including elitism (whereby the n most fit members of the population are automatically retained, and then the rest of the new generation is filled out randomly), and tournament selection (whereby two members of the population are picked out randomly, and then a further random number between 0 and 1 is selected.  If the random number exceeds a given parameter (often .25), the stronger member of the pair is chosen for the next generation; otherwise the weaker member is chosen.  Both members are then thrown back into the pool for potential re-selection.  The process is iterated for as many times as there are members of a generation).

But the classic method of selection, initially developed by John Holland, is "fitness-proportionate" selection.  Under fitness-proportionate selection, the odds that an individual will be selected to pass into the next generation are given as \(TraditionalForm
      \`\(individual' s\ fitness\)\/\(total\ fitness\ of\ population\)\).  There is no guarantee under this system that the strongest individual will move into the next generation, nor that the weakest individual will be extinguished.  Nevertheless, fitter individuals have a better than average chance to succeed.  Given the encoding we've chosen for our problem, fitness-proportionate selection will be implemented as follows.

First, the total fitness of the entire population is determined by summing the effective revenue values of each individual offer.  I believe totalfitness needs to remain defined as Set Delayed, because as population is updated with each generation, the total fitness will change (and, we hope, increase).  If totalfitness were defined as Set, that would work with the first generation, but would not allow for updating between generations.

population  =  Table[member,  {30}]

totalfitness := 
      Apply[Plus,  Map[effectiverevenue,  population]]

Then, a "roulette wheel" is created through a list of the cumulative sums of the fitness of the elements in the population.  This "roulette wheel" can be imagined to have a total circumference of precisely the total fitness of the population.  The 30 slots in the roulette wheel vary in width:  slot 1 is roulettewheel[[2]]-roulettewheel[[1]] units, which is the fitness of the first offer in the population.  Slot 2 is roulettewheel[[3]]-roulettewheel[[2]] units, which is the fitness of the second offer in the population.  If offer X generates twice as much effective revenue as offer Y (and thus is twice as fit), its slot on the roulette wheel will occupy twice as much space, and therefore will have twice as much chance of turning up when the wheel is spun.

roulettewheel := wheel  = 
        FoldList[Plus,  0, Map[effectiverevenue,  population]]

For as many times as there are offers in a generation (30 in the current implementation), the wheel is spun by selecting a random real number in the range from 0 to the total fitness of the population.  If, for example, the fifth time the wheel is spun, a number is generated which goes in the 9th slot of the roulette wheel (that is, it is greater than roulettewheel[[9]] but less than roulettewheel[[10]], then the ninth member of population becomes the fifth member of newgeneration.  Once newgeneration is completely populated, and the operations of crossover and mutation (to be discussed below) have been performed on it, population will be redefined as newgeneration so the process can start over again.

spinofthewheel  :=  selectionmark  =  
        Random[Real,  {0,  totalfitness}]

selection  :=  {roulettewheel;  
        For[m = 1,  
          m ≤ 30,  m++, {spinofthewheel; 
            For[n  =  1,  n ≤ 30,  n++,  
              If[selectionmark  > wheel[[n]]  &&  
                  selectionmark ≤ wheel[[n + 1]],  
                newgeneration[[m]]  =  
                  population [[n]]]]}]}

D.  Genetic Operators

1.  Crossover

As Mitchell explains , "crossover exchanges subparts of two chromosomes, roughly mimicking biological recombination between two single-chromosome ('haploid') organisms."  (IGA, p.3)  Without some mechanism such as crossover, all that could happen is that the population would gradually be filled with exact replicas of the fittest members of the original generation.  Crossover, by combining elements from adjacent members of the population, yields new offers that may have higher fitness than anything in the original population.

There are several ways in which crossover can be implemented (discussed in IGA, pp. 175-177).  Of these perhaps the simplest is single-point crossover.  With single point crossover, there are two elements of randomness:  whether the crossover will happen at all, and if so, at what point.  Much of the work done in genetic algorithms, as Mitchell notes, have centered on crossover probabilities of between .5 and .8.  I will be using a rate of .7.  Once it has been determined that crossover will occur, the next question is where.  If there are two offers, or organisms, each with five elements (let's call them AAAAA and BBBBB), then there are four internal crossover points.  If, for example, the second crossover point is used, then the two offers would swap elements after the second, resulting in AABBB and BBAAA.  If the fourth crossover point is used, the new offers will be AAAAB and BBBBA. 

The implementation for us is as follows:

newgeneration = Table[0, {30}];

crossovermechanism[x_]  :=  {prob = Random[Real];  
        If[prob  >  .3, { crosspoint = Random[Integer, {1, 4}];  
            population[[x]] = 
              Join[Take[newgeneration[[x]],  crosspoint],  
                Take[newgeneration[[x + 1]],  {crosspoint + 1,  5}]];  
            population[[x + 1]] = 
              Join[Take[newgeneration[[x + 1]],  crosspoint],  
                Take[newgeneration[[x]],  {crosspoint + 1, 
                    5}]]},  {population[[x]] = newgeneration[[x]],  
            population[[x + 1]] = 
              newgeneration[[x + 1]]}]}

crossover  :=  
      Flatten[Map[crossovermechanism[#] &,  Range[1, 30, 2]], 2]

D2.  Mutation

Point mutation is another operation which enables the "gene pool" to be somewhat transformed, thus creating the opportunity for organisms of higher fitness to arise.  Traditionally, rates of mutation are much lower than crossover rates (mimicking the lower incidence of point mutation in biological evolution).  We will use a fairly standard probability of 1 in 1000.  Since there are five different elements in each offer, each of which has somewhat different parameters, each needs a slightly different form of mutation.

For carrier and access method, mutation is simply a matter of changing the sign of the utility from positive to negative or vice versa

carriermutation[x_List] :=  {mutprob  =  Random[Real];  
        If[mutprob <  .001, MapAt[-# &,  x,  1],  
          MapAt[# &,  x,  1]]}

accessmethodmutation[x_List]  :=  {mutprob  =  Random[Real];  
        If[mutprob <  .001, MapAt[-# &,  x,  2],  
          MapAt[# &,  x,  2]]}

For the monthly fee, mutation will change the amount by a maximum of 25 cents.

monthlyfeemutation[x_List] :=  {mutprob  =  Random[Real]; 
        mutamount = Random[Real,  {- .25,   .25}];   If[
          mutprob  <   .001, MapAt[#  +  mutamount &,  x, 3],  
          MapAt[# &, x, 3]]}

For the per-call surcharge, mutation will change the amount by a maximum of 10 cents.

percallmutation[x_List] :=  {mutprob  =  Random[Real]; 
        mutamount = Random[Real,  {- .10,   .10}];   If[
          mutprob  <   .001, MapAt[#  +  mutamount &,  x, 4],  
          MapAt[# &, x, 4]]}

For the per-minute rate, mutation will change the amount by a maximum of 5 cents.

perminutemutation[x_List] :=  {mutprob  =  Random[Real]; 
        mutamount = Random[Real,  {- .05,   .05}];   If[
          mutprob  <   .001, MapAt[#  +  mutamount &,  x, 5],  
          MapAt[# &, x, 5]]}

mutation  :=  {population  = 
          Flatten[Map[carriermutation,  population], 1];  
        population  =  Flatten[Map[accessmethodmutation,  population], 1]; 
        population  = Flatten[ Map[monthlyfeemutation,  population], 1];  
        population  = Flatten[Map[percallmutation,  population], 1]; 
        population  = 
          Flatten[ Map[perminutemutation,  population], 1]}

E.  Putting it all together

All of the key elements and processes have now been defined.  The flow of evolution is as follows:

1.  We begin with a randomly generated population of thirty offers (population), a list of 30 initially empty elements into which can be mapped the members of the next generation (newgeneration), a generation counter set to 1 (generation), and lists of  initially empty elements into which can be mapped the value of the highest fitness achieved by any offer within a generation (maxfitness) and the average fitness of all offers within a generation (avgfitness).
2.  The selection process unfolds as described in section IVc above, with the result that newgeneration is created from the operation of the selection process upon population.
3.  The crossover process unfolds as described in section IVd1, with the result that a revised population is created from the operation of crossover on newgeneration.
4.  The mutation process unfolds as described in section IVd2, potentially modifying population.
5.  The generation counter is incremented.
6.  The fitness value of the offer yielding the greatest effective revenue from population is mapped into maxfitness[[generation]], while the average fitness of all offers in population is mapped onto avgfitness[[generation]]
7.  Steps 2 through 6 are repeated any desired number of times.  

Thus:

population  =  Table[member,  {30}]

We benchmark the results by computing the maximum fitness within, and the average fitness of, the randomly given starting population.

startingmaxfitness  = 
      effectiverevenue[
        Part[Sort[population,  
            effectiverevenue[#1] > effectiverevenue[#2] &], 1]]

startingavgfitness  =  totalfitness/30

newgeneration =  Table[0, {30}]; maxfitness  =  Table[0,  {50}];  
    avgfitness =  Table[0, {50}];

evolution  :=  {selection,  crossover,  mutation}

Note:  On my Pentium 166 computer, this takes between 1 and 2 minutes to run

For[generation  = 1,  
      generation  < 51,  generation++,  {evolution; 
        maxfitness[[generation]] = 
          effectiverevenue[
            Part[Sort[population,  
                effectiverevenue[#1] > effectiverevenue[#2] &], 1]];  
        avgfitness[[generation]] =  totalfitness/30}];  
    Show[{ListPlot[
          avgfitness,  {PlotStyle -> {Red},  PlotJoined -> True}],  
        ListPlot[
          maxfitness,  {PlotStyle -> {Blue},  
            PlotJoined -> True}]}]

Kevin's original results

bestevolvedoffer  = 
      Part[Sort[population,  effectiverevenue[#1] > effectiverevenue[#2] &], 
        1]

effectiverevenue[bestevolvedoffer]

Conclusions

The above run is fairly typical.  It turns out that, with the given parameters, the genetic algorithm does only a modest job of improving the maximum fitness within the population (usually only about 10%-15%above the maximum fitness within the randomly given population).  Moreover, most of the improvement in the maximum fitness takes place within the first few generations.  Because fitter options are more likely to be selected to reproduce, however, the average fitness of the population usually improves considerably, very nearly approaching the maximum fitness within 20 generations or so. 

It will be worth exploring in the futurewhether the performance of the algorithm in increasing maximum fitness can be improved by such means as increasing the mutation rate, changing the selection process to prevent early convergence of solutions, or adopting parameterized single-bit crossover  (see Mitchell, IGA)