diff --git a/README.md b/README.md
index f6b317e..856acc9 100644
--- a/README.md
+++ b/README.md
@@ -17,7 +17,7 @@ All of the algorithms are parallelized with Rayon.
 You need to implement your own model that implements `OptModel` trait. Actual optimization is handled by each algorithm functions. Here is a simple example to optimize a quadratic function with Hill Climbing algorithm.
 
 ```rust
-use std::error::Error;
+use std::{error::Error, time::Duration};
 
 use indicatif::{ProgressBar, ProgressDrawTarget, ProgressStyle};
 use localsearch::{
@@ -98,6 +98,7 @@ fn main() {
 
     println!("running Hill Climbing optimizer");
     let n_iter = 10000;
+    let time_limit = Duration::from_secs(60);
     let patiance = 1000;
     let n_trials = 50;
     let opt = HillClimbingOptimizer::new(patiance, n_trials);
@@ -107,10 +108,11 @@ fn main() {
         pb.set_position(op.iter as u64);
     };
 
-    let res = opt.optimize(&model, None, n_iter, Some(&callback), ());
+    let res = opt.optimize(&model, None, n_iter, time_limit, Some(&callback), ());
     pb.finish();
     dbg!(res);
 }
+
 ```
 
 Further details can be found at API document, example and test codes.
\ No newline at end of file
diff --git a/examples/quadratic_model.rs b/examples/quadratic_model.rs
index 6496c0e..86dd574 100644
--- a/examples/quadratic_model.rs
+++ b/examples/quadratic_model.rs
@@ -1,4 +1,4 @@
-use std::error::Error;
+use std::{error::Error, time::Duration};
 
 use indicatif::{ProgressBar, ProgressDrawTarget, ProgressStyle};
 use localsearch::{
@@ -79,6 +79,7 @@ fn main() {
 
     println!("running Hill Climbing optimizer");
     let n_iter = 10000;
+    let time_limit = Duration::from_secs_f32(1.0);
     let patiance = 1000;
     let n_trials = 50;
     let opt = HillClimbingOptimizer::new(patiance, n_trials);
@@ -88,7 +89,7 @@ fn main() {
         pb.set_position(op.iter as u64);
     };
 
-    let res = opt.optimize(&model, None, n_iter, Some(&callback), ());
+    let res = opt.optimize(&model, None, n_iter, time_limit, Some(&callback), ());
     pb.finish();
     dbg!(res);
 }
diff --git a/examples/tsp_model.rs b/examples/tsp_model.rs
index 47876f2..740226a 100644
--- a/examples/tsp_model.rs
+++ b/examples/tsp_model.rs
@@ -1,3 +1,4 @@
+use core::time::Duration;
 use std::{
     collections::{HashMap, HashSet},
     error::Error,
@@ -243,6 +244,7 @@ fn main() {
     let tsp_model = TSPModel::from_coords(&coords);
 
     let n_iter: usize = 20000;
+    let time_limit = Duration::from_secs(60);
     let patience = n_iter / 2;
 
     let mut rng = rand::thread_rng();
@@ -266,6 +268,7 @@ fn main() {
         &tsp_model,
         initial_solution.clone(),
         n_iter,
+        time_limit,
         Some(&callback),
         (),
     );
@@ -284,6 +287,7 @@ fn main() {
         &tsp_model,
         initial_solution.clone(),
         n_iter,
+        time_limit,
         Some(&callback),
         tabu_list,
     );
@@ -301,6 +305,7 @@ fn main() {
         &tsp_model,
         initial_solution.clone(),
         n_iter,
+        time_limit,
         Some(&callback),
         (200.0, 50.0),
     );
@@ -318,6 +323,7 @@ fn main() {
         &tsp_model,
         initial_solution.clone(),
         n_iter,
+        time_limit,
         Some(&callback),
         (),
     );
@@ -331,8 +337,14 @@ fn main() {
 
     println!("run relative annealing");
     let optimizer = RelativeAnnealingOptimizer::new(patience, 200, 10, 1e1);
-    let (final_solution, final_score, _) =
-        optimizer.optimize(&tsp_model, initial_solution, n_iter, Some(&callback), ());
+    let (final_solution, final_score, _) = optimizer.optimize(
+        &tsp_model,
+        initial_solution,
+        n_iter,
+        time_limit,
+        Some(&callback),
+        (),
+    );
     println!(
         "final score = {}, num of cities {}",
         final_score,
diff --git a/src/optim/base.rs b/src/optim/base.rs
index 6098a17..4ebc048 100644
--- a/src/optim/base.rs
+++ b/src/optim/base.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use auto_impl::auto_impl;
 use trait_set::trait_set;
 
@@ -17,6 +19,7 @@ pub trait LocalSearchOptimizer<M: OptModel> {
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
diff --git a/src/optim/epsilon_greedy.rs b/src/optim/epsilon_greedy.rs
index f84053e..e3f98eb 100644
--- a/src/optim/epsilon_greedy.rs
+++ b/src/optim/epsilon_greedy.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use crate::{callback::OptCallbackFn, OptModel};
 
 use super::{base::LocalSearchOptimizer, GenericLocalSearchOptimizer};
@@ -53,6 +55,7 @@ impl<M: OptModel> LocalSearchOptimizer<M> for EpsilonGreedyOptimizer {
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         _extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
@@ -66,6 +69,13 @@ impl<M: OptModel> LocalSearchOptimizer<M> for EpsilonGreedyOptimizer {
             self.return_iter,
             |current, trial| transition_prob(current, trial, self.epsilon),
         );
-        optimizer.optimize(model, initial_solution, n_iter, callback, _extra_in)
+        optimizer.optimize(
+            model,
+            initial_solution,
+            n_iter,
+            time_limit,
+            callback,
+            _extra_in,
+        )
     }
 }
diff --git a/src/optim/generic.rs b/src/optim/generic.rs
index cb031a3..54d944f 100644
--- a/src/optim/generic.rs
+++ b/src/optim/generic.rs
@@ -1,4 +1,9 @@
-use std::{cell::RefCell, marker::PhantomData, rc::Rc};
+use std::{
+    cell::RefCell,
+    marker::PhantomData,
+    rc::Rc,
+    time::{Duration, Instant},
+};
 
 use rand::Rng;
 use rayon::prelude::*;
@@ -69,12 +74,14 @@ where
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         _extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
     where
         F: OptCallbackFn<M::SolutionType, M::ScoreType>,
     {
+        let start_time = Instant::now();
         let mut rng = rand::thread_rng();
         let mut current_solution = if let Some(s) = initial_solution {
             s
@@ -88,6 +95,10 @@ where
         let mut counter = 0;
 
         for it in 0..n_iter {
+            let duration = Instant::now().duration_since(start_time);
+            if duration > time_limit {
+                break;
+            }
             let (trial_solution, trial_score) = (0..self.n_trials)
                 .into_par_iter()
                 .map(|_| {
diff --git a/src/optim/hill_climbing.rs b/src/optim/hill_climbing.rs
index 7ae9b42..4279777 100644
--- a/src/optim/hill_climbing.rs
+++ b/src/optim/hill_climbing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use crate::{callback::OptCallbackFn, OptModel};
 
 use super::{EpsilonGreedyOptimizer, LocalSearchOptimizer};
@@ -34,6 +36,7 @@ impl<M: OptModel> LocalSearchOptimizer<M> for HillClimbingOptimizer {
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         _extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
@@ -41,6 +44,13 @@ impl<M: OptModel> LocalSearchOptimizer<M> for HillClimbingOptimizer {
         F: OptCallbackFn<M::SolutionType, M::ScoreType>,
     {
         let optimizer = EpsilonGreedyOptimizer::new(self.patience, self.n_trials, usize::MAX, 0.0);
-        optimizer.optimize(model, initial_solution, n_iter, callback, _extra_in)
+        optimizer.optimize(
+            model,
+            initial_solution,
+            n_iter,
+            time_limit,
+            callback,
+            _extra_in,
+        )
     }
 }
diff --git a/src/optim/logistic_annealing.rs b/src/optim/logistic_annealing.rs
index 01ccae8..59472d6 100644
--- a/src/optim/logistic_annealing.rs
+++ b/src/optim/logistic_annealing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use ordered_float::NotNan;
 
 use crate::{callback::OptCallbackFn, OptModel};
@@ -58,6 +60,7 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for LogisticA
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         _extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
@@ -71,7 +74,14 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for LogisticA
             |current, trial| transition_prob(current, trial, self.w),
         );
 
-        optimizer.optimize(model, initial_solution, n_iter, callback, _extra_in)
+        optimizer.optimize(
+            model,
+            initial_solution,
+            n_iter,
+            time_limit,
+            callback,
+            _extra_in,
+        )
     }
 }
 
diff --git a/src/optim/relative_annealing.rs b/src/optim/relative_annealing.rs
index 2773e03..003644b 100644
--- a/src/optim/relative_annealing.rs
+++ b/src/optim/relative_annealing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use ordered_float::NotNan;
 
 use crate::{callback::OptCallbackFn, OptModel};
@@ -59,6 +61,7 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for RelativeA
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         _extra_in: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
@@ -72,7 +75,14 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for RelativeA
             |current, trial| transition_prob(current, trial, self.w),
         );
 
-        optimizer.optimize(model, initial_solution, n_iter, callback, _extra_in)
+        optimizer.optimize(
+            model,
+            initial_solution,
+            n_iter,
+            time_limit,
+            callback,
+            _extra_in,
+        )
     }
 }
 
diff --git a/src/optim/simulated_annealing.rs b/src/optim/simulated_annealing.rs
index b04cb13..d043a28 100644
--- a/src/optim/simulated_annealing.rs
+++ b/src/optim/simulated_annealing.rs
@@ -1,4 +1,8 @@
-use std::{cell::RefCell, rc::Rc};
+use std::{
+    cell::RefCell,
+    rc::Rc,
+    time::{Duration, Instant},
+};
 
 use ordered_float::NotNan;
 use rand::Rng;
@@ -46,12 +50,14 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for Simulated
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         max_min_temperatures: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
     where
         F: OptCallbackFn<M::SolutionType, M::ScoreType>,
     {
+        let start_time = Instant::now();
         let (max_temperature, min_temperature) = max_min_temperatures;
         let mut rng = rand::thread_rng();
         let mut current_solution = if let Some(s) = initial_solution {
@@ -68,6 +74,10 @@ impl<M: OptModel<ScoreType = NotNan<f64>>> LocalSearchOptimizer<M> for Simulated
         let mut counter = 0;
 
         for it in 0..n_iter {
+            let duration = Instant::now().duration_since(start_time);
+            if duration > time_limit {
+                break;
+            }
             let (trial_solution, trial_score) = (0..self.n_trials)
                 .into_par_iter()
                 .map(|_| {
diff --git a/src/optim/tabu_search.rs b/src/optim/tabu_search.rs
index 9879150..441bc82 100644
--- a/src/optim/tabu_search.rs
+++ b/src/optim/tabu_search.rs
@@ -1,4 +1,9 @@
-use std::{cell::RefCell, marker::PhantomData, rc::Rc};
+use std::{
+    cell::RefCell,
+    marker::PhantomData,
+    rc::Rc,
+    time::{Duration, Instant},
+};
 
 use auto_impl::auto_impl;
 use rayon::prelude::*;
@@ -91,12 +96,14 @@ impl<M: OptModel, T: TabuList<Item = (M::SolutionType, M::TransitionType)>> Loca
         model: &M,
         initial_solution: Option<M::SolutionType>,
         n_iter: usize,
+        time_limit: Duration,
         callback: Option<&F>,
         mut tabu_list: Self::ExtraIn,
     ) -> (M::SolutionType, M::ScoreType, Self::ExtraOut)
     where
         F: OptCallbackFn<M::SolutionType, M::ScoreType>,
     {
+        let start_time = Instant::now();
         let mut rng = rand::thread_rng();
         let mut current_solution = if let Some(s) = initial_solution {
             s
@@ -110,6 +117,10 @@ impl<M: OptModel, T: TabuList<Item = (M::SolutionType, M::TransitionType)>> Loca
         let mut accepted_counter = 0;
 
         for it in 0..n_iter {
+            let duration = Instant::now().duration_since(start_time);
+            if duration > time_limit {
+                break;
+            }
             let mut samples = vec![];
             (0..self.n_trials)
                 .into_par_iter()
diff --git a/src/tests/test_epsilon_greedy.rs b/src/tests/test_epsilon_greedy.rs
index 18c07c1..db7b397 100644
--- a/src/tests/test_epsilon_greedy.rs
+++ b/src/tests/test_epsilon_greedy.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::optim::{EpsilonGreedyOptimizer, LocalSearchOptimizer};
@@ -9,7 +11,14 @@ fn test() {
     let model = QuadraticModel::new(3, vec![2.0, 0.0, -3.5], (-10.0, 10.0));
     let opt = EpsilonGreedyOptimizer::new(1000, 10, 200, 0.1);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) = opt.optimize(&model, None, 10000, null_closure, ());
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        10000,
+        Duration::from_secs(10),
+        null_closure,
+        (),
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.05);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.05);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.05);
diff --git a/src/tests/test_hill_climbing.rs b/src/tests/test_hill_climbing.rs
index b95c31e..37202bf 100644
--- a/src/tests/test_hill_climbing.rs
+++ b/src/tests/test_hill_climbing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::optim::{HillClimbingOptimizer, LocalSearchOptimizer};
@@ -9,7 +11,14 @@ fn test() {
     let model = QuadraticModel::new(3, vec![2.0, 0.0, -3.5], (-10.0, 10.0));
     let opt = HillClimbingOptimizer::new(1000, 10);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) = opt.optimize(&model, None, 10000, null_closure, ());
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        10000,
+        Duration::from_secs(10),
+        null_closure,
+        (),
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.05);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.05);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.05);
diff --git a/src/tests/test_logistic_annealing.rs b/src/tests/test_logistic_annealing.rs
index a3685d0..3810165 100644
--- a/src/tests/test_logistic_annealing.rs
+++ b/src/tests/test_logistic_annealing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::optim::{LocalSearchOptimizer, LogisticAnnealingOptimizer};
@@ -9,7 +11,14 @@ fn test() {
     let model = QuadraticModel::new(3, vec![2.0, 0.0, -3.5], (-10.0, 10.0));
     let opt = LogisticAnnealingOptimizer::new(5000, 10, 200, 1e1);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) = opt.optimize(&model, None, 10000, null_closure, ());
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        10000,
+        Duration::from_secs(10),
+        null_closure,
+        (),
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.05);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.05);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.05);
diff --git a/src/tests/test_relative_annealing.rs b/src/tests/test_relative_annealing.rs
index d78b766..5245aae 100644
--- a/src/tests/test_relative_annealing.rs
+++ b/src/tests/test_relative_annealing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::optim::{LocalSearchOptimizer, RelativeAnnealingOptimizer};
@@ -9,7 +11,14 @@ fn test() {
     let model = QuadraticModel::new(3, vec![2.0, 0.0, -3.5], (-10.0, 10.0));
     let opt = RelativeAnnealingOptimizer::new(5000, 10, 200, 1e1);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) = opt.optimize(&model, None, 10000, null_closure, ());
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        10000,
+        Duration::from_secs(10),
+        null_closure,
+        (),
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.05);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.05);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.05);
diff --git a/src/tests/test_simulated_annealing.rs b/src/tests/test_simulated_annealing.rs
index ea6df77..5526dee 100644
--- a/src/tests/test_simulated_annealing.rs
+++ b/src/tests/test_simulated_annealing.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::optim::{LocalSearchOptimizer, SimulatedAnnealingOptimizer};
@@ -9,8 +11,14 @@ fn test() {
     let model = QuadraticModel::new(3, vec![2.0, 0.0, -3.5], (-10.0, 10.0));
     let opt = SimulatedAnnealingOptimizer::new(10000, 10);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) =
-        opt.optimize(&model, None, 5000, null_closure, (1.0, 0.1));
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        5000,
+        Duration::from_secs(10),
+        null_closure,
+        (1.0, 0.1),
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.05);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.05);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.05);
diff --git a/src/tests/test_tabu_search.rs b/src/tests/test_tabu_search.rs
index 1b07f44..4c5c801 100644
--- a/src/tests/test_tabu_search.rs
+++ b/src/tests/test_tabu_search.rs
@@ -1,3 +1,5 @@
+use std::time::Duration;
+
 use approx::assert_abs_diff_eq;
 
 use crate::{
@@ -40,8 +42,14 @@ fn test() {
     let opt = TabuSearchOptimizer::new(1000, 25, 5);
     let tabu_list = MyTabuList::new(10);
     let null_closure = None::<&fn(_)>;
-    let (final_solution, final_score, _) =
-        opt.optimize(&model, None, 10000, null_closure, tabu_list);
+    let (final_solution, final_score, _) = opt.optimize(
+        &model,
+        None,
+        10000,
+        Duration::from_secs(10),
+        null_closure,
+        tabu_list,
+    );
     assert_abs_diff_eq!(2.0, final_solution[0], epsilon = 0.1);
     assert_abs_diff_eq!(0.0, final_solution[1], epsilon = 0.1);
     assert_abs_diff_eq!(-3.5, final_solution[2], epsilon = 0.1);