diff --git a/mealpy/__init__.py b/mealpy/__init__.py
index b5ba005f..941f7669 100644
--- a/mealpy/__init__.py
+++ b/mealpy/__init__.py
@@ -128,3 +128,4 @@ def get_optimizer_by_name(name, verbose=False):
         for name, optimizer in cls.items():
             print(f"Optimizer: {name} - {optimizer} - {optimizer()}")
     return cls
+from .RFO import OriginalRFO
diff --git a/mealpy/swarm_based/RFO.py b/mealpy/swarm_based/RFO.py
new file mode 100644
index 00000000..81481deb
--- /dev/null
+++ b/mealpy/swarm_based/RFO.py
@@ -0,0 +1,33 @@
+import numpy as np
+from mealpy.optimizer import Optimizer
+
+
+class OriginalRFO(Optimizer):
+    """
+    Red Fox Optimization (RFO)
+    Dehghani et al., Artificial Intelligence Review, 2021
+    DOI: 10.1007/s10462-020-09904-6
+    """
+
+    def __init__(self, epoch=1000, pop_size=50, **kwargs):
+        super().__init__(**kwargs)
+        self.epoch = epoch
+        self.pop_size = pop_size
+        self.set_parameters(["epoch", "pop_size"])
+        self.sort_flag = False
+
+    def evolve(self, epoch):
+        a = 2 * (1 - epoch / self.epoch)
+        pop_new = []
+
+        for i in range(self.pop_size):
+            r = np.random.uniform(-1, 1, self.problem.n_dims)
+            dist = np.abs(self.g_best.solution - self.pop[i].solution)
+            pos_new = self.g_best.solution + a * r * dist
+            pos_new = self.amend_solution(pos_new)
+
+            agent = self.generate_empty_agent(pos_new)
+            agent.target = self.get_target(pos_new)
+            pop_new.append(agent)
+
+        self.pop = self.update_target_for_population(pop_new)
diff --git a/mealpy/test_RFO.py b/mealpy/test_RFO.py
new file mode 100644
index 00000000..1ecce8a9
--- /dev/null
+++ b/mealpy/test_RFO.py
@@ -0,0 +1,21 @@
+import numpy as np
+from mealpy import FloatVar
+from mealpy.swarm_based.RFO import OriginalRFO
+
+
+def sphere(x):
+    return np.sum(x ** 2)
+
+
+def test_rfo_runs():
+    problem = {
+        "obj_func": sphere,
+        "bounds": FloatVar(lb=(-5.,) * 5, ub=(5.,) * 5),
+        "minmax": "min",
+    }
+
+    model = OriginalRFO(epoch=10, pop_size=10)
+    result = model.solve(problem)
+
+    assert result.solution is not None
+
diff --git a/tests/swarm_based/test_RFO.py b/tests/swarm_based/test_RFO.py
new file mode 100644
index 00000000..73e2184c
--- /dev/null
+++ b/tests/swarm_based/test_RFO.py
@@ -0,0 +1,21 @@
+import numpy as np
+from mealpy import FloatVar
+from mealpy.swarm_based import OriginalRFO
+
+
+def objective_function(solution):
+    return np.sum(solution ** 2)
+
+
+def test_OriginalRFO():
+    problem = {
+        "obj_func": objective_function,
+        "bounds": FloatVar(lb=(-10.,) * 5, ub=(10.,) * 5),
+        "minmax": "min",
+    }
+
+    model = OriginalRFO(epoch=50, pop_size=20)
+    result = model.solve(problem)
+
+    assert result.solution is not None
+    assert result.target.fitness is not None