I'm using [Grammatical Evolution](https://en.wikipedia.org/wiki/Grammatical_evolution) for [symbolic regression](https://en.wikipedia.org/wiki/Symbolic_regression) tasks (i.e. searching the space of mathematical expressions to find the model that best fits a given dataset).

I'm evolving solutions according to a user-specified grammar:

```
 (0)
 <expr> ::= <expr> <op> <expr>      (0)
          | <const>                 (1)
          | <var>                   (2)
          | <unop>                  (3)

 (1)     
 <op> ::= *     (0)
        | +     (1)
        | -     (2) 
        | /     (3)

 (2)
 <const> ::= 0.1            (0)
           | 1.0            (1)
           | 10.0           (2)

 (3)
 <var> ::= x                (0) 

 (4)
 <unop> ::= sin             (0)
          | cos             (1) 
          | -               (2)

```

The population is a collection of arrays of integers. For example:

```
[0, 1, 2, 0, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1]

```

is the genotype of the individual `+(x, +(1, +(1, 1))` (function $x + 3$).

The translation is performed using the integers to look up the table which describe the [BNF](https://en.wikipedia.org/wiki/Backus%E2%80%93Naur_Form) grammar:

```
Step  Sequence  Rule   State
  1      0       0.0   <expr> <op> <expr>
  2      1       1.1   <expr> + <expr>
  3      2       0.2   <var> + <expr>
  4      0       3.0   x + <expr>
  5      0       0.0   x + <expr> <op> <expr>
  6      1       1.1   x + <expr> + <expr>
  7      1       0.1   x + <const> + <expr>
  8      1       0.1   x + 1.0 + <expr>
  9      0       0.1   x + 1.0 + <expr> <op> <expr>
 10      1       1.1   x + 1.0 + <expr> + <expr>
 11      1       0.1   x + 1.0 + <const> + <expr>
 12      1       2.1   x + 1.0 + 1.0 + <expr>
 13      1       0.1   x + 1.0 + 1.0 + <const>
 14      1       2.1   x + 1.0 + 1.0 + 1.0

```

The fitness function is the reciprocal of the sum, taken over the fitness
cases, of the absolute error between the evolved and the target function.

My problem is that, while the algorithm finds the general model, it often has trouble finding the right constants.

For example, if the target function is $x^2 + x - 12$, it finds something like:

$$x^2 + x - (10 + \sin(x))$$

or

$$x^2 + x - (1 + 1) \times (10 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1) \times ((10 - (10 / (1 + 1)) - (10 / (1 + 1 + 1 + 1 + 1)))$$

which is correct but requires a great effort to discover the constant.

Are there efficient techniques for generating numerical constants?