diff --git a/Project.toml b/Project.toml
index 75614acf..99b7ebd9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Metaheuristics"
 uuid = "bcdb8e00-2c21-11e9-3065-2b553b22f898"
 authors = ["Jesus Mejia <jesusmejded@gmail.com>"]
-version = "3.2.14"
+version = "3.2.15"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
diff --git a/docs/src/api.md b/docs/src/api.md
index a65644eb..5a0ef53b 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -106,6 +106,10 @@ Metaheuristics.compare
 Metaheuristics.gen_initial_state
 ```
 
+```@docs
+set_user_solutions!
+```
+
 ## Variation 
 
 ```@docs
diff --git a/src/DecisionMaking/ROI.jl b/src/DecisionMaking/ROI.jl
index dac61a41..92e524f6 100644
--- a/src/DecisionMaking/ROI.jl
+++ b/src/DecisionMaking/ROI.jl
@@ -91,7 +91,7 @@ function roiarchiving(population, w, parameters::ROIArchiving; verbose=true)
     
     δ_w = parameters.δ_w_transformed
 
-    size(w,2) != length(δ_w) && error("|weight_points| is different to |δ_w|")
+    size(w,1) != length(δ_w) && error("|weight_points| is different to |δ_w|")
 
     non_dominated = Metaheuristics.get_non_dominated_solutions_perm(population)
     isempty(non_dominated) && return Int[]
diff --git a/src/Metaheuristics.jl b/src/Metaheuristics.jl
index 267aad54..4da88a80 100644
--- a/src/Metaheuristics.jl
+++ b/src/Metaheuristics.jl
@@ -26,6 +26,7 @@ export GenerationalReplacement, ElitistReplacement
 export sample, LatinHypercubeSampling, Grid 
 export εDE, Restart
 export optimize!
+export set_user_solutions!
 
 include("externals.jl")
 
@@ -51,6 +52,8 @@ include("common/repair.jl")
 include("common/stop.jl")
 include("common/compare.jl")
 include("common/non-dominated-sorting.jl")
+include("common/set_user_solutions.jl")
+
 
 
 include("TestProblems/TestProblems.jl")
diff --git a/src/algorithms/ABC/ABC.jl b/src/algorithms/ABC/ABC.jl
index cd1c3676..4f7f4218 100644
--- a/src/algorithms/ABC/ABC.jl
+++ b/src/algorithms/ABC/ABC.jl
@@ -90,7 +90,10 @@ function initialize!(
         options.debug && @info "Increasing f calls limit to $(options.f_calls_limit)"
     end
 
-    bees = initialbees(parameters.N, problem)
+
+    _st = gen_initial_state(problem,parameters,information,options,status)
+    bees = [Bee(sol) for sol in _st.population]
+
     nevals = length(bees)
 
     best_sol = deepcopy(getBestBee(bees))
diff --git a/src/algorithms/ABC/bee_dynamics.jl b/src/algorithms/ABC/bee_dynamics.jl
index 9ee951b5..6489408b 100644
--- a/src/algorithms/ABC/bee_dynamics.jl
+++ b/src/algorithms/ABC/bee_dynamics.jl
@@ -137,9 +137,10 @@ function chooseBest(bees, best)
     return best
 end
 
+#=
 function initialbees(N, problem)
      P = generate_population(N, problem)
 
      return [ Bee(sol) for sol in P ]
 end 
-
+=#
diff --git a/src/algorithms/CGSA/CGSA.jl b/src/algorithms/CGSA/CGSA.jl
index 45f2be27..10a911e1 100644
--- a/src/algorithms/CGSA/CGSA.jl
+++ b/src/algorithms/CGSA/CGSA.jl
@@ -135,17 +135,11 @@ function initialize!(
 
     max_it = 500
     options.iterations = options.iterations == 0 ? max_it : options.iterations
-    options.f_calls_limit = options.f_calls_limit == 0 ? options.iterations * N : options.f_calls_limit
 
     # random initialization for agents.
-    P = generate_population(N, problem)
-
-    # Current best
-    theBest = get_best(P)
-
-
-    status = State(theBest, P)
-    status.f_calls = N
+    status = gen_initial_state(problem,parameters,information,options,status)
+    N = parameters.N
+    options.f_calls_limit = options.f_calls_limit == 0 ? options.iterations * N : options.f_calls_limit
 
     # Velocity
     parameters.V = isempty(parameters.V) ? zeros(N,D) : parameters.V
@@ -153,8 +147,7 @@ function initialize!(
     parameters.X = isempty(parameters.X) ? positions(status) : parameters.X
     # function values
     if isempty(parameters.fitness)
-
-        parameters.fitness = fval.(status.population)
+        parameters.fitness = fvals(status.population)
     end
     
 
diff --git a/src/algorithms/MCCGA/MCCGA.jl b/src/algorithms/MCCGA/MCCGA.jl
index e6d34043..0aafc47a 100644
--- a/src/algorithms/MCCGA/MCCGA.jl
+++ b/src/algorithms/MCCGA/MCCGA.jl
@@ -125,11 +125,8 @@ function initialize!(
 
     parameters.probvector = initialprobs(lower, upper, maxsamples = parameters.maxsamples)
 
-    # sample a vector to create an initial State
-    x = sample(parameters.probvector) |> floats
-    initial_sol = create_solution(x, problem)
-    return State(initial_sol, [initial_sol for i in 1:parameters.N])
-
+    # sample vectors to create an initial State
+    return gen_initial_state(problem,parameters,information,options,status)
 end
 
 function update_state!(
diff --git a/src/algorithms/WOA/WOA.jl b/src/algorithms/WOA/WOA.jl
index 66230098..3b826ef1 100644
--- a/src/algorithms/WOA/WOA.jl
+++ b/src/algorithms/WOA/WOA.jl
@@ -82,16 +82,14 @@ function initialize!(
   options.iterations * parameters.N : options.f_calls_limit
 
 
-  P = generate_population(parameters.N, problem)
-  best_sol = deepcopy(get_best(P))
-  status = State(best_sol, P)
+  # P = generate_population(parameters.N, problem)
+  # best_sol = deepcopy(get_best(P))
+  # status = State(best_sol, P)
   #=
   status.population = P
   status.best_sol = deepcopy(get_best(status.population))
   =#
-  status.f_calls = parameters.N
-
-  status
+  return gen_initial_state(problem,parameters,information,options,status)
 
 end
 
diff --git a/src/algorithms/algorithm.jl b/src/algorithms/algorithm.jl
index 824f2f10..e69de29b 100644
--- a/src/algorithms/algorithm.jl
+++ b/src/algorithms/algorithm.jl
@@ -1,42 +0,0 @@
-"""
-    gen_initial_state(problem,parameters,information,options)
-
-Generate an initial state, i.e., compute uniformly distributed random vectors in bounds,
-after that are evaluated in objective function. This method require that `parameters.N`
-is valid attribute.
-"""
-gen_initial_state(problem,parameters,information,options,status::State{Any}) = gen_initial_state(problem,parameters,information,options)
-
-
-function gen_initial_state(problem,parameters,information,options, status)
-    parameters.N != length(status.population) &&
-        error("Population size in provided State differs from that in parameters")
-
-
-    size(problem.bounds,2) != length(get_position(status.best_sol)) &&
-        error("Invalid population (dimension does not match with bounds)")
-
-    return State(status.best_sol, status.population)
-
-    
-end
-
-function gen_initial_state(problem,parameters,information,options)
-    # population array
-    population = generate_population(parameters.N, problem,ε=options.h_tol)
-
-    # best solution
-    best_solution = get_best(population)
-
-    return State(best_solution, population; f_calls = length(population), iteration=1)
-end
-
-function Base.show(io::IO, parameters::AbstractParameters)
-    s = typeof(parameters)
-
-    vals = string.(map(f -> getfield(parameters, f), fieldnames(s)))
-    str = string(s) * "(" * join(string.(fieldnames(s)) .* "=" .* vals, ", ") * ")"
-
-    print(io, str)
-end
-
diff --git a/src/common/gen_initial_state.jl b/src/common/gen_initial_state.jl
new file mode 100644
index 00000000..75f8ab20
--- /dev/null
+++ b/src/common/gen_initial_state.jl
@@ -0,0 +1,67 @@
+"""
+    gen_initial_state(problem,parameters,information,options)
+
+Generate an initial state, i.e., compute uniformly distributed random vectors in bounds,
+after that are evaluated in objective function. This method require that `parameters.N`
+is valid attribute.
+"""
+gen_initial_state(problem,parameters,information,options,status::State{Any}) = gen_initial_state(problem,parameters,information,options)
+
+
+function gen_initial_state(problem,parameters,information,options, status)
+    if parameters.N != length(status.population)
+        options.debug && @warn("Population size in provided State differs from that in parameters")
+    end
+
+    size(problem.bounds,2) != length(get_position(status.best_sol)) &&
+        error("Invalid population (dimension does not match with bounds)")
+
+    _complete_population!(status,problem,parameters,information,options)
+
+    best_solution = get_best(status.population)
+    # check if a better solution was found
+    if is_better(status.best_sol, best_solution)
+        best_solution = status.best_sol
+    end
+
+    return State(best_solution, status.population;f_calls = length(status.population))
+end
+
+function gen_initial_state(problem,parameters,information,options)
+    # population array
+    population = generate_population(parameters.N, problem,ε=options.h_tol)
+
+    # best solution
+    best_solution = get_best(population)
+
+    return State(best_solution, population; f_calls = length(population), iteration=1)
+end
+
+function Base.show(io::IO, parameters::AbstractParameters)
+    s = typeof(parameters)
+
+    vals = string.(map(f -> getfield(parameters, f), fieldnames(s)))
+    str = string(s) * "(" * join(string.(fieldnames(s)) .* "=" .* vals, ", ") * ")"
+
+    print(io, str)
+end
+
+function _complete_population!(status,problem,parameters,information,options)
+    if parameters.N < length(status.population)
+        # increase population if necessary
+        parameters.N = length(status.population)
+        # TODO: use this options.debug == true to show the message?
+        @warn("Population size increased to $(parameters.N) due to initial solutions.")
+        return
+    end
+
+    if parameters.N > length(status.population)
+        # complete population with random in bounds
+        n = parameters.N - length(status.population)
+        missing_sols = generate_population(n, problem,ε=options.h_tol)
+        # insert new solution into population
+        append!(status.population, missing_sols)
+    end
+    
+end
+
diff --git a/src/common/set_user_solutions.jl b/src/common/set_user_solutions.jl
new file mode 100644
index 00000000..757d78b1
--- /dev/null
+++ b/src/common/set_user_solutions.jl
@@ -0,0 +1,86 @@
+include("gen_initial_state.jl")
+
+"""
+    set_user_solutions!(optimizer, x, fx)
+
+Provide initial solutions to the `optimizer`.
+- `x` can be a `Vector` and `fx` a function or `fx = f(x)`
+- `x` can be a matrix containing solutions in rows.
+
+### Example
+
+```julia-repl
+julia> f(x) = abs(x[1]) + x[2]  + x[3]^2 # objective function
+f (generic function with 1 method)
+
+julia> algo  = ECA(N = 61); # optimizer
+
+julia> # one solution can be provided
+       x0 = [0.5, 0.5, 0.5];
+
+julia> set_user_solutions!(algo, x0, f);
+
+julia> # providing multiple solutions
+       X0 = rand(30, 3); # 30 solutions with dim 3
+
+julia> set_user_solutions!(algo, X0, f);
+
+julia> optimize(f, [0 0 0; 1 1 1.0], algo)
++=========== RESULT ==========+
+  iteration: 413
+    minimum: 0
+  minimizer: [0.0, 0.0, 0.0]
+    f calls: 25132
+ total time: 0.0856 s
+stop reason: Small difference of objective function values.
++============================+
+```
+
+"""
+function set_user_solutions!(algo::AbstractAlgorithm, solution::AbstractSolution)
+    status = algo.status
+    if !isnothing(status.best_sol) && !isempty(status.population)
+        push!(status.population, solution)
+    else
+        algo.status = State(solution, [solution])
+    end
+    algo
+end
+
+
+function set_user_solutions!(algo::AbstractAlgorithm, x::AbstractVector, fx)
+    set_user_solutions!(algo, create_child(x, fx))
+end
+
+function set_user_solutions!(algo::AbstractAlgorithm, x::AbstractVector, f::Function)
+    set_user_solutions!(algo, x, f(x))
+end
+
+
+function set_user_solutions!(algo::AbstractAlgorithm, X::AbstractMatrix, fX::AbstractVector)
+    n = size(X, 1)
+    m = length(fX)
+
+    if n != m
+        @warn "$(n) decision vectors provided but $(m) objective values."
+        n = min(m, n)
+        println("Taking ", n, " as the number of initial solutions.")
+    end
+
+    # nothing to do due to it is necessary the objective value
+    n == 0 && (return algo)
+
+    # TODO: this part can be parallelized
+    for i in 1:n
+        set_user_solutions!(algo, X[i,:], fX[i])
+    end
+
+    # TODO check population size provided in algo.parameters.N
+
+    algo
+end
+
+function set_user_solutions!(algo::AbstractAlgorithm, X::AbstractMatrix, f::Function)
+    set_user_solutions!(algo, X, [f(X[i,:]) for i in 1:size(X,1)])
+end
+
diff --git a/test/decisionmaking.jl b/test/decisionmaking.jl
index 2a903005..ed3f059a 100644
--- a/test/decisionmaking.jl
+++ b/test/decisionmaking.jl
@@ -21,20 +21,20 @@
 end
 
 @testset "DecisionMaking: ROI" begin
-    function test_roi()
+    function test_roi(w = [0.1 0.9; 0.9 0.1], δ = [0.1,0.1])
         _, _, pf = Metaheuristics.TestProblems.ZDT1();
         res = State(pf[1], pf)
-
-        w = [0.1 0.9; 0.9 0.1]
-        δ = [0.1,0.1]
+        
         method = ROIArchiving(δ)
         idx = decisionmaking(res, w, method)
         @test !isempty(idx)
         sols = best_alternative(res, w, method)
         @test length(sols) == length(idx)
     end
+
     
     test_roi() 
+    test_roi([0.0 1; 1 0; 0.5 0.5], [0.1, 0.2, 0.3]) 
 end
 
 @testset "DecisionMaking: JMcDM" begin
diff --git a/test/initial_solution.jl b/test/initial_solution.jl
new file mode 100644
index 00000000..6ceab1ed
--- /dev/null
+++ b/test/initial_solution.jl
@@ -0,0 +1,80 @@
+@testset "Initial Solutions" begin
+
+    function run_initialized(problem)
+        ff, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)
+
+        # number of function evaluations
+        f_calls = 0
+
+        f(x) = begin
+            f_calls += 1
+            ff(x)
+        end
+
+        options = Options(iterations = 3, f_calls_limit = 500)
+        information = Information()
+        x_optimum = get_position(first(optimums))
+
+        methods = [ABC, DE, ECA, PSO, WOA, MCCGA, CGSA]
+
+        for method in methods
+            for N in [5, 10, 20] # vary the population size
+                algo = method(N = N, options = options, information = information)
+                # initial solutions
+                X = rand(10, size(bounds, 2))
+                X[1,:] = x_optimum #zeros(size(bounds, 2))
+                set_user_solutions!(algo, X, f)
+
+                res = optimize(f, bounds, algo)
+
+                @test length(res.population) == algo.parameters.N 
+                @test minimum(res) ≈ first(fvals(optimums))
+                @test sum(abs,minimizer(res)-x_optimum) ≈ 0
+                # TODO check the number of function evaluations
+                #@test f_calls == nfes(res)
+            end
+        end
+    end
+
+    function run_initialized_multiobjective(problem)
+        ff, bounds, optimums = Metaheuristics.TestProblems.get_problem(problem)
+
+        # number of function evaluations
+        f_calls = 0
+
+        f(x) = begin
+            f_calls += 1
+            ff(x)
+        end
+
+        options = Options(iterations = 3, f_calls_limit = 500)
+        information = Information()
+        x_optimum = get_position(first(optimums))
+
+        methods = [SPEA2, SMS_EMOA, NSGA2, NSGA3]
+
+        for method in methods
+            for N in [5, 10, 20] # vary the population size
+                algo = method(N = N, options = options, information = information)
+                # initial solutions
+                X = rand(10, size(bounds, 2))
+                X[1,:] = x_optimum
+                set_user_solutions!(algo, X, f)
+
+                res = optimize(f, bounds, algo)
+                @test length(res.population) == algo.parameters.N 
+                # check whether optimum was used by optimizer
+                v = sum(abs.(positions(res) .- x_optimum'), dims=2) |> minimum 
+                @test  v ≈ 0
+            end
+        end
+    end
+
+    for problem in [:sphere]
+        run_initialized(problem)
+    end
+
+    for problem in [:ZDT1]
+        run_initialized_multiobjective(problem)
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 51c2d505..7e312ce4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,7 +16,8 @@ for tests in [
               "constrained.jl",
               "multi-objective.jl",
               "combinatorial.jl",
-              "decisionmaking.jl"
+              "decisionmaking.jl",
+              "initial_solution.jl",
              ]
     include(tests)
 end