Genetic algorithm does not converge to exact solution

Question

I'm trying to evolve quantum circuits using genetic algorithms as they did in this paper Decomposition of unitary matrices for finding quantum circuits: Application to molecular Hamiltonians (Daskin & Kais).

So given a unitary target matrix ($U_t$), the task is to evolve a circuit that approximates $U_t$ (let's call it $U_a$). This is the function I'm maximizing:

$$F = \frac{1}{2^n}|\operatorname{Tr}(U_aU_t^{\dagger})|$$ (where n is the number of qubits)

My circuit encoding is a little bit different, I'm using a json string representation, here's an example of a random circuit:

[{0: ('Y', -1, None), 1: ('Td', -1, None), 2: ('Y', -1, None)}, 
{0: ('Z', 2, None), 1: ('Y', -1, None)}, 
{0: ('Vd', 1, None), 2: ('V', -1, None)}, 
{1: ('X', -1, None), 2: ('Y', 0, None)}, 
{0: ('Td', -1, None), 2: ('V', 1, None)}, 
{1: ('S', -1, None), 2: ('X', 0, None)}]

Every dictionary item in the array is a column in the format {target: (gate, control, param)}. The value -1 means no control, and None means there are no parameters. Here's a visual representation of this circuit:

+--------+-------+-------+-------+-------+-------+-------+
| qubit  | col_0 | col_1 | col_2 | col_3 | col_4 | col_5 |
+--------+-------+-------+-------+-------+-------+-------+
| q0-->  |   Y   |  Z•2  |  Vd•1 |   I   |   Td  |   I   |
| q1-->  |   Td  |   Y   |   I   |   X   |   I   |   S   |
| q2-->  |   Y   |   I   |   V   |  Y•0  |  V•1  |  X•0  |
+--------+-------+-------+-------+-------+-------+-------+

So far, I've been successful at finding very small circuits (with 2 "columns"). The problem comes when trying to evolve solutions for more complex circuits. For example, the paper shows this solution for the Toffoli unitary:

My algorithm is not able to find this solution or any other exact solution to this problem. I tried changing the selection method, crossover and mutation rates, but the fitness value never exceeds 0.85. After some generations, it seems that all individuals are converging to the same non-optimal circuit (this happens even with uniform crossover and random selection).

For crossover, I've tried one-point crossover and uniform crossover. The mutation operation can do these things:

1. If it's a controlled gate, it can change the type of the target gate, can change the position of the control, can swap target/control, can remove the control, or can shuffle a column.
2. If it's not a controlled gate, it can change the type of gate or introduce a control.

I don't see what else can I do to improve these operators. In the paper, they also have a continuous angle phase parameter, but this particular circuit solution does not use it anyway.

Another curious fact: for this particular problem, the empty circuit has a fitness of 0.75. If a take the exact solution (from the image) and change a single gate, I can get a fitness of 0.5. So, an empty circuit has better fitness than a circuit which is a small step from the exact solution. This is very confusing, I don't see where the GA will find "optimal sub-circuits" to exchange during crossover.

It should be noted that the paper uses a more sophisticated version of a GA, but for the reasons stated above, I don't think it will make a difference in my case. So what I'm missing or how can I improve this? Thank you!

Edit

It seems the GLOA did it's work. The algorithm found this for the Toffoli unitary:

+--------+-------+-------+-------+-------+-------+
| qubit  | col_0 | col_1 | col_2 | col_3 | col_4 |
+--------+-------+-------+-------+-------+-------+
| q0-->  |   I   |   I   |   I   |   I   |   I   |
| q1-->  |  Z•2  |   I   |   I   |   I   |  S•0  |
| q2-->  |   I   |  Vd•0 |  Z•1  |  V•0  |   I   |
+--------+-------+-------+-------+-------+-------+

But I limited the kind of gates it could use, now I'll try with all the gates.