diff --git a/core/MarkovProcesses.jl b/core/MarkovProcesses.jl
index 4ae10bbef973bcd676c37a17fe4201c6ef979429..946d155f36bb24e0928ee995ff53710196bcaef9 100644
--- a/core/MarkovProcesses.jl
+++ b/core/MarkovProcesses.jl
@@ -26,7 +26,7 @@ export isbounded, isaccepted, check_consistency
 export load_model, get_module_path
 
 # Utils
-export get_module_path
+export get_module_path, cosmos_get_values
 
 include("common.jl")
 
diff --git a/core/utils.jl b/core/utils.jl
index 4128525ddd046e04f908a87b609048dd7f9ce5a3..5b46b3bcd1d53d10f1302766a4be967c6293dcda 100644
--- a/core/utils.jl
+++ b/core/utils.jl
@@ -1,3 +1,17 @@
 
 get_module_path() = dirname(dirname(pathof(@__MODULE__)))
 
+function cosmos_get_values(name_file::String) 
+    output_file = open(name_file)
+    dict_values = Dict{String}{Float64}()
+    while !eof(output_file)
+        line = readline(output_file)
+        splitted_line = split(line, ':')
+        if (length(splitted_line) > 1) && tryparse(Float64, splitted_line[2]) !== nothing 
+            dict_values[splitted_line[1]] = parse(Float64, splitted_line[2])
+        end
+    end
+    close(output_file)
+    return dict_values
+end
+
diff --git a/models/ER.jl b/models/ER.jl
index deace15ceccebd55fbdfb4d550b3ca867ce15cc4..3c4a3c467b483d3f8159bf2a79f55e58e77fbf16 100644
--- a/models/ER.jl
+++ b/models/ER.jl
@@ -46,5 +46,13 @@ isabsorbing_ER(p::Vector{Float64},xn::SubArray{Int,1}) =
 g_ER = ["P"]
 
 ER = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_ER,p_ER,x0_ER,t0_ER,ER_f!,isabsorbing_ER; g=g_ER)
-export ER
+
+function create_ER(new_p::Vector{Float64})
+    ER_new = deepcopy(ER)
+    @assert length(ER_new.p) == length(new_p)
+    set_param!(ER_new, new_p)
+    return ER_new
+end
+
+export ER, create_ER
 
diff --git a/models/SIR.jl b/models/SIR.jl
index 37990bca6de67b84490b84e01964f55a4c23882f..c9033d45aa7e5888bbeb3668dd8aea430310a741 100644
--- a/models/SIR.jl
+++ b/models/SIR.jl
@@ -45,5 +45,13 @@ isabsorbing_SIR(p::Vector{Float64}, xn::SubArray{Int,1}) = (p[1]*xn[1]*xn[2] + p
 g_SIR = ["I"]
 
 SIR = ContinuousTimeModel(d,k,dict_var,dict_p,l_tr_SIR,p_SIR,x0_SIR,t0_SIR,SIR_f!,isabsorbing_SIR; g=g_SIR)
-export SIR
+
+function create_SIR(new_p::Vector{Float64})
+    SIR_new = deepcopy(SIR)
+    @assert length(SIR_new.p) == length(new_p)
+    set_param!(SIR_new, new_p)
+    return SIR_new
+end
+
+export SIR, create_SIR
 
diff --git a/tests/cosmos/distance_F/ER.gspn b/tests/cosmos/distance_F/ER.gspn
new file mode 100644
index 0000000000000000000000000000000000000000..bd6462eed78cb5bd771fdec8468b3e84c58fbf07
--- /dev/null
+++ b/tests/cosmos/distance_F/ER.gspn
@@ -0,0 +1,45 @@
+const int E_init=100;
+const int S_init=100;
+const int ES_init=0;
+const int P_init=0;
+
+const double k_1=1.0;
+const double k_2=1.0;
+const double k_3=1.0;
+
+
+NbPlaces = 4;
+NbTransitions = 3;
+
+PlacesList = { 
+E,S, ES, P
+} ;
+
+TransitionsList = { 
+TR1,   TR2,  TR3
+} ;
+
+Marking={
+(E , E_init); (S , S_init); (ES , ES_init);
+(P , P_init); 
+};
+
+Transitions={
+(TR1,EXPONENTIAL(k_1*E*S),1,1, SINGLE);
+(TR2,EXPONENTIAL(k_2*ES),1,1, SINGLE);
+(TR3,EXPONENTIAL(k_3*ES),1,1, SINGLE);
+};
+
+InArcs={
+(E, TR1,1); (S, TR1,1); 
+(ES, TR2,1); (ES, TR3,1); 
+};
+
+OutArcs={
+(TR1, ES,1); 
+(TR2, E,1);
+(TR2, S,1); 
+(TR3, E,1); 
+(TR3, P,1); 
+};
+
diff --git a/tests/cosmos/distance_F/ER_1D.jl b/tests/cosmos/distance_F/ER_1D.jl
new file mode 100644
index 0000000000000000000000000000000000000000..4c7a50288c2f4106e3497ed7a336644a949e7d4c
--- /dev/null
+++ b/tests/cosmos/distance_F/ER_1D.jl
@@ -0,0 +1,71 @@
+
+
+using MarkovProcesses
+absolute_path = get_module_path() * "/tests/cosmos/"
+
+# Remarque: la valeur estimée par Cosmos varie plus que de 0.01
+
+# Values x1, x2  t1, t2
+dict_exp = Dict(
+                "R1" => [50, 75, 0.025, 0.05],
+                "R2" => [50, 75, 0.05, 0.075],
+                "R3" => [25, 50, 0.05, 0.075]
+               )
+
+str_model = "ER"
+if str_model == "ER"
+    load_model(str_model)
+    model = ER
+end
+load_automaton("automaton_F")
+
+test_all = true
+width = 0.1
+level = 0.95
+l_exp = ["R1", "R2", "R3"]
+#l_exp = ["R2"]
+for exp in l_exp
+    if !(exp in keys(dict_exp))
+        println("Unrecognized experiment: <<$exp>>")
+        continue
+    end
+    val_exp = dict_exp[exp]
+    x1, x2, t1, t2 = val_exp[1], val_exp[2], val_exp[3], val_exp[4]
+    A_F = create_automaton_F(model, x1, x2, t1, t2, "P")  
+    l_k3 = 0:10:100
+    nb_param = length(l_k3)
+    l_dist_cosmos = zeros(nb_param)
+    l_dist_pkg = zeros(nb_param)
+    for i in 1:nb_param
+        # Cosmos estimation
+        k3 = l_k3[i]
+        command = `Cosmos $(absolute_path * "distance_F/" * str_model * ".gspn") 
+        $(absolute_path * "distance_F/dist_F_"  * str_model * ".lha") --njob $(ENV["JULIA_NUM_THREADS"]) 
+        --const k_3=$(k3),x1=$x1,x2=$x2,t1=$t1,t2=$t2 
+        --level $(level) --width $(width) 
+        --verbose 0` 
+        run(pipeline(command, stderr=devnull))
+        dict_values = cosmos_get_values("Result_dist_F_$(str_model).res")
+        l_dist_cosmos[i] = dict_values["Estimated value"]
+        nb_sim = dict_values["Total paths"]
+        nb_accepted = dict_values["Accepted paths"]
+        nb_sim = convert(Int, nb_sim)
+        # MarkovProcesses estimation
+        set_param!(ER, "k3", convert(Float64, k3))
+        sync_ER = ER*A_F
+        l_dist_i_threads = zeros(Threads.nthreads()) 
+        Threads.@threads for j = 1:nb_sim
+            σ = simulate(sync_ER)
+            l_dist_i_threads[Threads.threadid()] += (σ.S)["d"]
+        end
+        l_dist_pkg[i] = sum(l_dist_i_threads) / nb_sim
+        global test_all = test_all && isapprox(l_dist_cosmos[i], l_dist_pkg[i]; atol = width)
+    end
+
+end
+
+rm("Result_dist_F_$(str_model).res")
+rm("Result.res")
+
+return test_all
+
diff --git a/tests/cosmos/distance_F/dist_F_ER.lha b/tests/cosmos/distance_F/dist_F_ER.lha
new file mode 100644
index 0000000000000000000000000000000000000000..baa9f54eb7d36ed1364653a087e4bddb77a3f050
--- /dev/null
+++ b/tests/cosmos/distance_F/dist_F_ER.lha
@@ -0,0 +1,40 @@
+NbVariables = 3;
+NbLocations = 4;
+
+const int x1 = 50;
+const int x2 = 75;
+
+const double t1 = 0.025;
+const double t2 = 0.05;
+
+VariablesList = {t,d,n};
+
+LocationsList = {l0, l1, l2, l3};
+
+AVG(Last(d));
+
+InitialLocations={l0};
+FinalLocations={l2};
+
+Locations={
+(l0,   TRUE , (t:1));
+(l1,   TRUE , (t:1));
+(l2,   TRUE , (t:1));
+(l3,   TRUE , (t:1));
+};
+
+Edges={
+((l0,l1), #, t>=0, {n=P});
+((l1,l3), #, t<=t1 & n<=x1-1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
+((l1,l3), #, t<=t1 & n>=x2+1, {d=min(((t-t1)^2+(n-x2)^2)^0.5,((t-t1)^2+(n-x1)^2)^0.5)});
+((l1,l3), #, n>=x1 & n<=x2, {d=0});
+((l1,l3), #, n<=x1-1 & t >= t1 & t<=t2, {d=min(d ,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) )});
+((l1,l3), #, n>=x2+1 & t >= t1 & t<=t2, {d=min(d,min( ((n-x1)^2)^0.5, ((n-x2)^2)^0.5 ) ) });
+((l3,l1), ALL, t>=0, {n=P});
+
+((l3,l2), #, t>=t2 ,#);
+
+((l1,l2),#, n>=x1 & n<=x2 & t>=t1 & t<=t2 ,{d=0});
+((l1,l2),#, t>=t2 , #);
+((l1,l2),#, d=0 & t>=t1, #);
+};
diff --git a/tests/run_all.jl b/tests/run_all.jl
index 5ed6f292c64e5b11d3a0008f10150b3ec9ec2cf0..67ed81726a22b2e5c547240f98702da7bcd84052 100644
--- a/tests/run_all.jl
+++ b/tests/run_all.jl
@@ -3,4 +3,5 @@ include("run_unit.jl")
 include("run_simulation.jl")
 include("run_dist_lp.jl")
 include("run_automata.jl")
+include("run_cosmos.jl")
 
diff --git a/tests/run_cosmos.jl b/tests/run_cosmos.jl
new file mode 100644
index 0000000000000000000000000000000000000000..096a7f53738bc03bcea540fe5426810b88b3a6c5
--- /dev/null
+++ b/tests/run_cosmos.jl
@@ -0,0 +1,7 @@
+
+using Test
+
+@testset "Cosmos tests" begin
+    @test include("cosmos/distance_F/ER_1D.jl")
+end
+
diff --git a/tests/run_unit.jl b/tests/run_unit.jl
index 49dfe1f240d1901d6d3945d4d99d216fe3d8e575..ea30b82044f1d2db764825e58ac47460a0e430c0 100644
--- a/tests/run_unit.jl
+++ b/tests/run_unit.jl
@@ -20,5 +20,6 @@ using Test
     @test include("unit/check_model_consistency.jl")
     @test include("unit/set_param.jl")
     @test include("unit/side_effects_1.jl")
+    @test include("unit/create_models.jl")
 end
 
diff --git a/tests/unit/create_models.jl b/tests/unit/create_models.jl
new file mode 100644
index 0000000000000000000000000000000000000000..1d2df65596f4be6086b75c620406c9ab5cd96ee8
--- /dev/null
+++ b/tests/unit/create_models.jl
@@ -0,0 +1,26 @@
+
+using MarkovProcesses
+
+load_model("SIR")
+load_model("ER")
+
+test_all = true
+
+new_p = [0.01, 0.2]
+new_SIR = create_SIR(new_p)
+test_all = test_all && new_SIR.p == new_p && new_SIR !== SIR
+
+p_SIR = [0.0012, 0.05]
+new_SIR = create_SIR(p_SIR)
+test_all = test_all && new_SIR.p == p_SIR && new_SIR !== SIR
+
+new_p = [0.01, 0.2, 20.0]
+new_ER = create_ER(new_p)
+test_all = test_all && new_ER.p == new_p && new_ER !== ER
+
+p_ER = [1.0, 1.0, 1.0]
+new_ER = create_ER(p_ER)
+test_all = test_all && new_ER.p == p_ER && new_ER !== ER
+
+return true
+