Evolutionary programming (EP) [Fogel92][Baeck93][Baeck97] is a computational technique that mimics evolution and is based on reproduction and selection. An EP algorithm is composed of individuals that reproduce and compete, each one is a potential solution to the (optimization) problem and is represented by a "genome" where each gene corresponds to one adjustable parameter. At each generation of the EP, each individual reproduces asexually, i.e. divides into two individuals. One of these contains exactly the same "genome" as the parent while the other suffers some mutations (the parameter values of each gene change slightly). At the end of the generation, the algorithm has double the number of individuals. Then each of the individuals is confronted with a number of others to count how many does it outperform (the number of wins is the number of these competitors that represent worse solutions than itself). All the individuals are ranked by their number of wins, and the population is again reduced to the original number of individuals by eliminating those which have worse fitness (solutions).

Evolutionary Strategies with Stochastic Ranking (SRES) [Runarsson00] is similar to Evolutionary Programming. However, a parent has multiple offsprings during each generation. Each offspring will contain a recombination of genes with another parent and additional mutations. The algorithm assures that each parameter value will be within its boundaries. But constraints to the solutions may be violated. Whenever this happens the square of the size of the violation is summed up, i.e., we calculate

where the constraints are given by . The value is used within the selection, which is performed as described in Genetic Algorithm SR

The genetic algorithm (GA) [Baeck97][Baeck93][Michalewicz94][Mitchell95] is a computational technique that mimics evolution and is based on reproduction and selection. A GA is composed of individuals that reproduce and compete, each one is a potential solution to the (optimization) problem and is represented by a "genome" where each gene corresponds to one adjustable parameter. At each generation of the GA, each individual is paired with one other at random for reproduction. Two offspring are produced by combining their genomes and allowing for "cross-over", i.e., the two new individuals have genomes that are formed from a combination of the genomes of their parents. Also each new gene might have mutated, i.e. the parameter value might have changed slightly. At the end of the generation, the algorithm has double the number of individuals. Then each of the individuals is confronted with a number of others to count how many does it outperform (the number of wins is the number of these competitors that represent worse solutions than itself). All the individuals are ranked by their number of wins, and the population is again reduced to the original number of individuals by eliminating those which have worse fitness (solutions).

Many features of a GA may be varied. The details of this particular implementation of the GA for optimization of biochemical kinetics are:

The genetic algorithm with stochastic ranking is very similar to the before described Genetic Algorithm with tournament selection. With two exception which are the mutations are not forced to be within the boundaries and the selection is done through a bubble sort with a random factor as described in [Runarsson00].

The method of Hooke and Jeeves [Bell66][Hooke61][Kaupe63][Swann72] is a direct search algorithm that searches for the minimum of a nonlinear function without requiring (or attempting to calculate) derivatives of the function. Instead it is based on a heuristic that suggests a descent direction using the values of the function calculated in a number of previous iterations.

Levenberg-Marquardt [Levenberg44][Marquardt63] is a gradient descent method. It is a hybrid between the steepest descent and the Newton methods.

The Newton optimization method searches for the minimum of a nonlinear function by following descent directions determined from the function's first and second partial derivatives. The steepest descent method searches for a minimum based only on the first derivatives of the function. While the Newton method converges quadratically towards the minimum in its vicinity, it may not converge at all if it is far away from it. On the other hand the steepest descent method only converges linearly but is guaranteed to converge.

Levenberg first suggested an improvement to the Newton method in order to make it more robust, i.e. to overcome the problem of non-convergence. His suggestion was to add a factor to the diagonal elements of the Hessian matrix of second derivatives when not close to the minimum (this can be judged by how positive definite the matrix is). The effect when this factor is large compared to the elements of Hessian is that the method then becomes the steepest descent method. Later Marquardt suggested that the factor should be multiplicative rather than additive and also defined a heuristic to make this factor increase or decrease. The method known as Levenberg-Marquardt is thus an adaptive method that effectively changes between the steepest descent to the Newton method.

The original suggestions of Levenberg and Marquardt were effective to enhance the Gauss-Newton method, a variant of the Newton method specifically for minimizing least-squares functions. In this case the advantage is also that the second derivatives do not need to be calculated as they are estimated from the gradient of the residuals. Subsequently Goldfeld et al. [Goldfeld66] extended the method to the case of general non-linear functions.

This method also known as the simplex method is due to Nelder and Mead [Nelder65]. A simplex is a polytope of vertices in dimensions. The objective function is evaluated at each vertex. Dependent on these calculated values a new simplex is constructed. The simplest step is to replace the worst point with a point reflected through the centroid of the remaining points. If this point is better than the best current point, then we can try stretching exponentially out along this line. On the other hand, if this new point isn't much better than the previous value then we are stepping across a valley, so we shrink the simplex towards the best point.

The particle swarm optimization method suggested by Kennedy and Eberhart [Kennedy95] is inspired by a flock of birds or a school of fish searching for food. Each particle has a position and a velocity in the parameter space. Additionally, it remembers its best achieved objective value and position . Dependent on its own information and the position of its best neighbor (a random subset of particles of the swarm) a new velocity is calculated. With this information the position is updated. The pseudo code for the algorithm for each particle is:

  N =  best neighbor's position

  for i = 0 to number of parameters do
    R1 = uniform random number in [0, 1]
    R2 = uniform random number in [0, 1]

    V[i]= w * V[i]
          + C * R1 * (M[i]- X[i])
          + C * R2 * (N[i]- X[i])
    X[i] = X[i] + V[i]
  enddo

  current_O = evaluate objective function

  if current_O < O then do
    O = current_O
    M = X
  enddo

Praxis is a direct search method (for a review see [Swann72]) that searches for the minimum of a nonlinear function without requiring (or attempting to calculate) derivatives of that function. Praxis was developed by Brent [Brent72] after the method proposed by Powell [Powel64]. The inspiration for Praxis was the well-known method of minimising each adjustable parameter (direction) at a time - the principal axes method. In Praxis directions are chosen that do not coincide with the principal axes, in fact if the objective function is quadratic then these will be conjugate directions, assuring a fast convergence rate.

This implementation of the Praxis method was originally written in FORTRAN at the Stanford Linear Accelerator Center (dated 3/1/73). The original FORTRAN code is available from http://www.netlib.org/opt/praxis. The original code was translated automatically by the F2C program from AT&T Bell Labs (available at http://www.netlib.org/f2c/). The C code was then enhanced to suite COPASI's reentrant optimization method interface.

Random search is an optimization method that attempts to find the optimum by testing the objective function's value on a series of combinations of random values of the adjustable parameters. The random values are generated complying with any boundaries selected by the user, furthermore, any combinations of parameter values that do not fulfill constraints on the variables are excluded. This means that the method is capable of handling bounds on the adjustable parameters and fulfilling constraints.

For infinite number of iterations this method is guaranteed to find the global optimum of the objective function. In general one is interested in processing a very large number of iterations.

Simulated annealing is an optimization algorithm first proposed by Kirkpatrick et al. [Kirkpatrick83] and was inspired by statistical mechanics and the way in which perfect crystals are formed. Perfect crystals are formed by first melting the substance of interest, and then cooling it very slowly. At large temperatures the particles vibrate with wide amplitude and this allows a search for global optimum. As the temperature decreases so do the vibrations until the system settles to the global optimum (the perfect crystal).

The simulated annealing optimization algorithm uses a similar concept: the objective function is considered a measure of the energy of the system and this is maintained constant for a certain number of iterations (a temperature cycle). In each iteration, the parameters are changed to a nearby location in parameter space and the new objective function value calculated; if it decreased, then the new state is accepted, if it increased then the new state is accepted with a probability that follows a Boltzmann distribution (higher temperature means higher probability of accepting the new state). After a fixed number of iterations, the stopping criterion is checked; if it is not time to stop, then the system's temperature is reduced and the algorithm continues.

Simulated annealing is a stochastic algorithm that is guaranteed to converge if ran for an infinite number of iterations. It is one of the most robust global optimization algorithms, although it is also one of the slowest. (Be warned that simulated annealing can run for hours or even days!).

This implementation of simulated annealing is based on the code of Corana et al. [Corana87]. This implementation has tuned some of the parameters of the original implementation as follows:

Steepest descent [Fogel92] is an optimization method that follows the direction of steepest descent on the hyper-surface of the objective function to find a local minimum. The direction of steepest descent is defined by the negative of the gradient of the objective function.

The Truncated Newton method is a sophisticated variant of the Newton optimisation method. The Newton optimisation method searches for the minimum of a nonlinear function by following descent directions determined from the function’s first and second partial derivatives. The Truncated Newton method does an incomplete (truncated) solution of a system of linear equations to calculate the Newton direction. This means that the actual direction chosen for the descent is between the steepest descent direction and the true Newton direction. A more complete description of the method can be found in [Gill81] and [Nash84].

The particular implementation of the Truncated Newton method used here takes into account simple bounds on the optimisation parameters. Like the Newton method itself, the Truncated Newton variant converges quadratically when in the vicinity of the minimum. On average this method is one of the best gradient descent methods.

This implementation of the Truncated Newton method was originally written in FORTRAN by Stephen Nash (Operations Research and Applied Statistics Dept., George Mason University, Fairfax, VA 22030, USA). The original code is available from http://www.netlib.org/opt/tn). The original code was translated automatically to C by the F2C program from AT&T Bell Labs (available at http://www.netlib.org/f2c/ ). The C code was then enhanced to suite COPASI's reentrant optimization method interface.