diff --git a/doc/guide/guide-control.rst b/doc/guide/guide-control.rst
index 915b6b6..09caedb 100644
--- a/doc/guide/guide-control.rst
+++ b/doc/guide/guide-control.rst
@@ -310,6 +310,107 @@ Eventually, the optimization for a desired `fid_err_targ` can be started by call
     },
   )
 
+The Genetic Algorithm (GA)
+==========================
+
+The Genetic Algorithm (GA) is a global optimization technique inspired by natural selection. 
+Unlike gradient-based methods like GRAPE or CRAB, GA explores the solution space stochastically, 
+making it robust against local minima and suitable for problems with non-differentiable or noisy objectives.
+
+In QuTiP, the GA-based optimizer evolves a population of candidate control pulses across multiple 
+generations to minimize the infidelity between the system's final and target states.
+
+The GA optimization consists of the following components:
+
+**Population**:
+    A collection of candidate solutions (chromosomes), where each chromosome encodes a full set 
+    of control amplitudes over time for the given control Hamiltonians.
+
+**Fitness Function**:
+    The fitness of each candidate is evaluated using a fidelity-based measure, such as:
+
+    - **PSU (Projective State Overlap)**: 
+      :math:`1 - |\langle \psi_{\text{target}} | \psi_{\text{final}} \rangle|`
+    - **TRACEDIFF**: 
+      Trace distance between final and target density matrices.
+
+**Genetic Operations**:
+
+- **Selection**:
+  A subset of the best-performing candidates (based on fitness) are chosen to survive.
+  
+- **Crossover (Mating)**:
+  New candidates are generated by combining genes from selected parents using arithmetic crossover.
+  
+- **Mutation**:
+  Random perturbations are added to the new candidates' genes to maintain diversity and escape local minima.
+
+This process continues until either a target fidelity error is reached or a fixed number of generations 
+have passed without improvement (stagnation criterion).
+
+Each generation represents a full evaluation of the population, making the method inherently parallelizable 
+and effective in high-dimensional control landscapes.
+
+In QuTiP, the GA optimization is implemented via the ``_Genetic`` class, and can be invoked using the 
+standard ``optimize_pulses`` interface by setting the algorithm to ``"Genetic"``.
+
+Optimal Quantum Control in QuTiP (Genetic Algorithm)
+====================================================
+
+Defining a control problem in QuTiP using the Genetic Algorithm follows the same pattern as with other algorithms.
+
+.. code-block:: python
+
+    import qutip as qt
+    import qutip_qoc as qoc
+    import numpy as np
+
+    # state to state transfer
+    initial = qt.basis(2, 0)
+    target = qt.basis(2, 1)
+
+    # define the drift and control Hamiltonians
+    drift = qt.sigmaz()
+    controls = [qt.sigmax(), qt.sigmay()]
+
+    H = [drift] + controls
+
+    # define the objective
+    objective = qoc.Objective(initial, H, target)
+
+    # discretized time grid (e.g., 100 steps over 10 units of time)
+    tlist = np.linspace(0, 10, 100)
+
+    # define control parameter bounds (same for all controls in GA)
+    control_parameters = {
+        "p": {"bounds": [(-1.0, 1.0)]}
+    }
+
+    # set genetic algorithm hyperparameters
+    algorithm_kwargs = {
+        "alg": "Genetic",
+        "population_size": 50,
+        "generations": 100,
+        "mutation_rate": 0.2,
+        "fid_err_targ": 1e-3,
+    }
+
+    # run the optimization
+    result = qoc.optimize_pulses(
+        objectives=[objective],
+        control_parameters=control_parameters,
+        tlist=tlist,
+        algorithm_kwargs=algorithm_kwargs,
+    )
+
+    print("Final infidelity:", result.infidelity)
+    print("Best control parameters:", result.optimized_params)
+
+References
+==========
+
+.. [Brown2023] J. Brown, M. Paternostro, and A. Ferraro, "Optimal quantum control via genetic algorithms for quantum state engineering in driven‑resonator mediated networks," *Quantum Sci. Technol.*, vol. 8, no. 2, p. 025004, 2023. https://doi.org/10.1088/2058-9565/acb2f2
+
 .. TODO: add examples
 
 Examples for Liouvillian dynamics and multi-objective optimization will follow soon.
diff --git a/src/qutip_qoc/_genetic.py b/src/qutip_qoc/_genetic.py
new file mode 100644
index 0000000..ede3ee8
--- /dev/null
+++ b/src/qutip_qoc/_genetic.py
@@ -0,0 +1,226 @@
+import numpy as np
+import qutip as qt
+import time
+from qutip_qoc.result import Result
+
+
+class _Genetic:
+    """
+    Genetic Algorithm-based optimizer for quantum control problems.
+
+    This class implements a global optimization routine using a Genetic Algorithm
+    to optimize control pulses that drive a quantum system from an initial state
+    to a target state (or unitary).
+
+    Based on the code from Jonathan Brown
+    """
+
+    def __init__(
+        self,
+        objectives,
+        control_parameters,
+        time_interval,
+        time_options,
+        alg_kwargs,
+        optimizer_kwargs,
+        minimizer_kwargs,
+        integrator_kwargs,
+        qtrl_optimizers,
+    ):
+        self._objective = objectives[0]
+        self._Hd = self._objective.H[0]
+        self._Hc_lst = [
+            H[0] if isinstance(H, list) else H for H in self._objective.H[1:]
+        ]
+        self._initial = self._objective.initial
+        self._target = self._objective.target
+        self._norm_fac = 1 / self._target.norm()
+
+        self._evo_time = time_interval.evo_time
+        self.N_steps = time_interval.n_tslots
+        self.N_controls = len(self._Hc_lst)
+        self.N_var = self.N_controls * self.N_steps
+
+        self._alg_kwargs = alg_kwargs
+        self.N_pop = alg_kwargs.get("population_size", 100)
+        self.generations = alg_kwargs.get("generations", 100)
+        self.mutation_rate = alg_kwargs.get("mutation_rate", 0.3)
+        self.fid_err_targ = alg_kwargs.get("fid_err_targ", 1e-4)
+        self._stagnation_patience = 50  # Internally fixed
+
+        self._integrator_kwargs = integrator_kwargs
+        self._fid_type = alg_kwargs.get("fid_type", "PSU")
+
+        self._generator = self._prepare_generator()
+        self._solver = (
+            qt.MESolver(H=self._generator, options=self._integrator_kwargs)
+            if self._Hd.issuper
+            else qt.SESolver(H=self._generator, options=self._integrator_kwargs)
+        )
+
+        self._result = Result(
+            objectives=[self._objective],
+            time_interval=time_interval,
+            start_local_time=time.time(),
+            guess_controls=[self._generator],
+            guess_params=control_parameters.get("guess"),
+            var_time=False,
+            qtrl_optimizers=qtrl_optimizers,
+        )
+        self._result.iter_seconds = []
+        self._result.infidelity = np.inf
+        self._result._final_states = []
+
+    def _prepare_generator(self):
+        args = {
+            f"p{i+1}_{j}": 0.0
+            for i in range(self.N_controls)
+            for j in range(self.N_steps)
+        }
+
+        def make_coeff(i, j):
+            return lambda t, args: (
+                args[f"p{i+1}_{j}"]
+                if int(t / (self._evo_time / self.N_steps)) == j
+                else 0
+            )
+
+        H_qev = [self._Hd]
+        for i, Hc in enumerate(self._Hc_lst):
+            for j in range(self.N_steps):
+                H_qev.append([Hc, make_coeff(i, j)])
+
+        return qt.QobjEvo(H_qev, args=args)
+
+    def _infid(self, params):
+        args = {
+            f"p{i+1}_{j}": params[i * self.N_steps + j]
+            for i in range(self.N_controls)
+            for j in range(self.N_steps)
+        }
+        result = self._solver.run(self._initial, [0.0, self._evo_time], args=args)
+        final_state = result.final_state
+        self._result._final_states.append(final_state)
+
+        if self._fid_type == "TRACEDIFF":
+            diff = final_state - self._target
+            fid = 0.5 * np.real((diff.dag() * diff).tr())
+        else:
+            overlap = self._norm_fac * self._target.overlap(final_state)
+            fid = (
+                1 - np.abs(overlap) if self._fid_type == "PSU" else 1 - np.real(overlap)
+            )
+
+        return fid
+
+    def step(self, params):
+        t0 = time.time()
+        val = -self._infid(params)
+        self._result.iter_seconds.append(time.time() - t0)
+        return val
+
+    def initial_population(self):
+        return np.random.uniform(-1, 1, (self.N_pop, self.N_var))
+
+    def darwin(self, population, fitness):
+        indices = np.argsort(-fitness)[: self.N_pop // 2]
+        return population[indices], fitness[indices]
+
+    def pairing(self, survivors, survivor_fitness):
+        prob_dist = survivor_fitness - np.min(survivor_fitness)
+        prob_dist /= np.sum(prob_dist)
+
+        mothers = np.random.choice(len(survivors), size=self.N_pop // 4, p=prob_dist)
+        fathers = []
+        for m in mothers:
+            p = prob_dist.copy()
+            p[m] = 0
+            p /= np.sum(p)
+            fathers.append(np.random.choice(len(survivors), p=p))
+        return survivors[mothers], survivors[fathers]
+
+    def mating_procedure(self, ma, da):
+        beta = np.random.rand(*ma.shape)
+        swap = np.random.randint(0, 2, ma.shape)
+        beta_inv = 1 - beta
+
+        new1 = swap * beta * ma + swap * beta_inv * da + (1 - swap) * ma
+        new2 = swap * beta * da + swap * beta_inv * ma + (1 - swap) * da
+        return np.vstack((new1, new2))
+
+    def build_next_gen(self, survivors, offspring):
+        return np.vstack((survivors, offspring))
+
+    def mutate(self, population):
+        n_mut = int(
+            (population.shape[0] - 1) * population.shape[1] * self.mutation_rate
+        )
+        row = np.random.randint(1, population.shape[0], size=n_mut)
+        col = np.random.randint(0, population.shape[1], size=n_mut)
+        population[row, col] += np.random.normal(0, 0.3, size=n_mut)
+        population[row, col] = np.clip(population[row, col], -2, 2)
+        return population
+
+    def optimize(self):
+        population = self.initial_population()
+        best_fit = -np.inf
+        best_chrom = None
+        history = []
+        no_improvement_counter = 0
+
+        for gen in range(self.generations):
+            fitness = np.array([self.step(chrom) for chrom in population])
+            max_fit = np.max(fitness)
+            history.append(max_fit)
+
+            if max_fit > best_fit:
+                best_fit = max_fit
+                best_chrom = population[np.argmax(fitness)]
+                no_improvement_counter = 0
+            else:
+                no_improvement_counter += 1
+
+            self._result.infidelity = min(self._result.infidelity, -max_fit)
+
+            if -max_fit <= self.fid_err_targ:
+                break
+
+            if no_improvement_counter >= self._stagnation_patience:
+                break
+
+            survivors, survivor_fit = self.darwin(population, fitness)
+            mothers, fathers = self.pairing(survivors, survivor_fit)
+            offspring = self.mating_procedure(mothers, fathers)
+            population = self.build_next_gen(survivors, offspring)
+            population = self.mutate(population)
+
+        self._result.optimized_params = best_chrom.tolist()
+        self._result.infidelity = -best_fit
+        self._result.end_local_time = time.time()
+        self._result.n_iters = len(history)
+        self._result.new_params = self._result.optimized_params
+        self._result._optimized_controls = best_chrom.tolist()
+        self._result._optimized_H = [self._generator]
+        self._result._final_states = self._result._final_states  # expose final_states
+        self._result.guess_params = self._result.guess_params or []
+        self._result.var_time = False
+
+        self._result.message = (
+            f"Stopped early: reached infidelity target {self.fid_err_targ}"
+            if -best_fit <= self.fid_err_targ
+            else (
+                f"Stopped due to stagnation after {self._stagnation_patience} generations"
+                if no_improvement_counter >= self._stagnation_patience
+                else "Optimization completed successfully"
+            )
+        )
+        return self._result
+
+    def result(self):
+        self._result.start_local_time = time.strftime(
+            "%Y-%m-%d %H:%M:%S", time.localtime(self._result.start_local_time)
+        )
+        self._result.end_local_time = time.strftime(
+            "%Y-%m-%d %H:%M:%S", time.localtime(self._result.end_local_time)
+        )
+        return self._result
diff --git a/src/qutip_qoc/pulse_optim.py b/src/qutip_qoc/pulse_optim.py
index 935bf18..6d7079f 100644
--- a/src/qutip_qoc/pulse_optim.py
+++ b/src/qutip_qoc/pulse_optim.py
@@ -12,6 +12,7 @@
 from qutip_qoc._time import _TimeInterval
 
 import qutip as qt
+from qutip_qoc._genetic import _Genetic
 
 try:
     from qutip_qoc._rl import _RL
@@ -418,6 +419,21 @@ def optimize_pulses(
         )
         rl_env.train()
         return rl_env.result()
+    
+    elif alg == "Genetic":
+        gen_env = _Genetic(
+            objectives,
+            control_parameters,
+            time_interval,
+            time_options,
+            algorithm_kwargs,
+            optimizer_kwargs,
+            minimizer_kwargs,
+            integrator_kwargs,
+            qtrl_optimizers,
+        )
+        gen_env.optimize()
+        return gen_env.result()
 
     return _global_local_optimization(
         objectives,
@@ -429,4 +445,4 @@ def optimize_pulses(
         minimizer_kwargs,
         integrator_kwargs,
         qtrl_optimizers,
-    )
+    )
\ No newline at end of file
diff --git a/tests/test_result.py b/tests/test_result.py
index 6bd7f8c..1043191 100644
--- a/tests/test_result.py
+++ b/tests/test_result.py
@@ -220,6 +220,41 @@ def sin_z_jax(t, r, **kwargs):
     state2state_rl = None
     unitary_rl = None
 
+# --------------------------- Genetic ---------------------------
+
+
+# state to state transfer
+init = qt.basis(2, 0)
+target = qt.basis(2, 1)
+
+H_c = [qt.sigmax(), qt.sigmay(), qt.sigmaz()] # control Hamiltonians
+
+w, d, y = 0.1, 1.0, 0.1
+H_d = 1 / 2 * (w * qt.sigmaz() + d * qt.sigmax()) # drift Hamiltonian
+
+H = [H_d] + H_c # total Hamiltonian
+
+state2state_genetic = Case(
+    objectives=[Objective(initial, H, target)],
+    control_parameters={
+        "p": {"bounds": [(-13, 13)]},
+    },
+    tlist=np.linspace(0, 10, 100), 
+    algorithm_kwargs={
+        "fid_err_targ": 0.01,
+        "alg": "Genetic",
+        "max_iter": 100,
+    },
+    optimizer_kwargs={},
+)
+
+
+initial = qt.qeye(2) # Identity
+target  = qt.gates.hadamard_transform()
+
+unitary_genetic = state2state_genetic._replace(
+    objectives=[Objective(initial, H, target)],
+)
 
 @pytest.fixture(
     params=[
@@ -229,6 +264,8 @@ def sin_z_jax(t, r, **kwargs):
         pytest.param(state2state_goat, id="State to state (GOAT)"),
         pytest.param(state2state_jax, id="State to state (JAX)"),
         pytest.param(state2state_rl, id="State to state (RL)"),
+        pytest.param(state2state_genetic, id="State to state (Genetic)"),
+        pytest.param(unitary_genetic, id="Unitary (Genetic)"),
         pytest.param(unitary_rl, id="Unitary (RL)"),
     ]
 )