diff --git a/BayesValidRox/surrogate_models/surrogate_models.py b/BayesValidRox/surrogate_models/surrogate_models.py
index e896dca3df166bcdca22f5a467063009c4f4d735..49bdb374b728b43c4ef6ca5ab64e37d4d7ca132f 100644
--- a/BayesValidRox/surrogate_models/surrogate_models.py
+++ b/BayesValidRox/surrogate_models/surrogate_models.py
@@ -1076,11 +1076,11 @@ class aPCE:
                 # find an optimal point subset to add to the initial design
                 # by maximization of the utility score and taking care of NaN values
                 temp = ExploitScore.copy()
-                bestIdx = np.argsort(temp)[::-1][:NrofNewSample]
-                #bestIdx = sorted_idxtotalScore[:NrofNewSample]
+                sorted_idxtotalScore = np.argsort(temp, axis=0)[::-1]
+                bestIdx = sorted_idxtotalScore[:NrofNewSample]
                 
                 Xnew = np.zeros((NrofNewSample,ndim))
-                for i, idx in enumerate(bestIdx):
+                for i, idx in enumerate(bestIdx.flatten()):
                     X_can = exploration.closestPoints[idx]
                     
                     #plotter(self.ExpDesign.X, X_can, ExploreMethod, scoreExploration=None)
diff --git a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
index 8143a81b9b73059b13ee3111e13d63c861c887ca..0628f19add72818f31e36e30e91d3edd6892169b 100644
--- a/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
+++ b/BayesValidRox/tests/AnalyticalFunction/AnalyticalFunction_Test.py
@@ -138,7 +138,7 @@ if __name__ == "__main__":
     #MetaModelOpts.ExpDesign.nReprications = 2
     # -------- Exploration ------
     #1)'Voronoi' 2)'MC' 3)'LHS' 4)'dual annealing'
-    MetaModelOpts.ExpDesign.ExploreMethod = 'MC'
+    MetaModelOpts.ExpDesign.ExploreMethod = 'Voronoi'
     
     # Use when 'dual annealing' chosen
     MetaModelOpts.ExpDesign.MaxFunItr = 200
@@ -149,16 +149,16 @@ if __name__ == "__main__":
     
     # -------- Exploitation ------
     # 1)'BayesOptDesign' 2)'VarOptDesign' 3)'alphabetic' 4)'Space-filling'
-    MetaModelOpts.ExpDesign.ExploitMethod = 'BayesOptDesign'
+    MetaModelOpts.ExpDesign.ExploitMethod = 'VarOptDesign'
     
     # BayesOptDesign -> when data is available
     # 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
     # 3)APP (A-Posterior-percision) 
-    MetaModelOpts.ExpDesign.UtilityFunction = 'DKL' #['DKL', 'DPP']
+    #MetaModelOpts.ExpDesign.UtilityFunction = 'DKL' #['DKL', 'DPP']
     
     # VarBasedOptDesign -> when data is not available
     # Only with Vornoi >>> 1)Entropy 2)EIGF, 3)ALM, 4)LOOCV
-    #MetaModelOpts.ExpDesign.UtilityFunction = 'Entropy'#['EIGF', 'Entropy', 'LOOCV']
+    MetaModelOpts.ExpDesign.UtilityFunction = 'LOOCV'#['EIGF', 'Entropy', 'LOOCV']
     
     # alphabetic
     # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)
diff --git a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
index e6959dd571dacb1274fbba7e6a5ec18ff37b42e6..4039df9c9ce92b7f3c2717859c20d23f41f7a904 100644
--- a/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
+++ b/BayesValidRox/tests/Co2BenchmarkTest/Test_CO2Benchmark.py
@@ -107,16 +107,16 @@ if __name__ == "__main__":
     MetaModelOpts.addExpDesign()
     MetaModelOpts.ExpDesign.MCSize = 10000
     
-    initNSamples = 200
+    initNSamples = 100
     MetaModelOpts.ExpDesign.NrSamples = initNSamples
-    MetaModelOpts.ExpDesign.SamplingMethod = 'user' # 1)LHS 2)PCM 3)LSCM 4)user
+    MetaModelOpts.ExpDesign.SamplingMethod = 'MC' # 1)LHS 2)PCM 3)LSCM 4)user
     MetaModelOpts.ExpDesign.Method = 'sequential'  # 1) normal  2) sequential
-    MetaModelOpts.ExpDesign.X = np.load('EDX_CO2Benchmark.npy')
+    #MetaModelOpts.ExpDesign.X = np.load('EDX_CO2Benchmark.npy')
     #MetaModelOpts.ExpDesign.Y = ExpDesign_Y
     
     # Sequential experimental design (needed only for sequential ExpDesign)
     MetaModelOpts.ExpDesign.NrofNewSample = 10
-    MetaModelOpts.ExpDesign.MaxNSamples = 500
+    MetaModelOpts.ExpDesign.MaxNSamples = 1000
     MetaModelOpts.ExpDesign.ModifiedLOOThreshold = 1e-4
     
     # Defining the measurement error, if it's known a priori
@@ -143,7 +143,7 @@ if __name__ == "__main__":
     #MetaModelOpts.ExpDesign.nReprications = 2
     # -------- Exploration ------
     #1)'Voronoi' 2)'MC' 3)'LHS' 4)'dual annealing'
-    MetaModelOpts.ExpDesign.ExploreMethod = 'MC'
+    MetaModelOpts.ExpDesign.ExploreMethod = 'Voronoi'
     
     # Use when 'dual annealing' chosen
     MetaModelOpts.ExpDesign.MaxFunItr = 200
@@ -154,16 +154,16 @@ if __name__ == "__main__":
     
     # -------- Exploitation ------
     # 1)'BayesOptDesign' 2)'VarOptDesign' 3)'alphabetic' 4)'Space-filling'
-    MetaModelOpts.ExpDesign.ExploitMethod = 'BayesOptDesign'
+    MetaModelOpts.ExpDesign.ExploitMethod = 'VarOptDesign'
     
     # BayesOptDesign -> when data is available
     # 1)DKL (Kullback-Leibler Divergence) 2)DPP (D-Posterior-percision)
     # 3)APP (A-Posterior-percision) 
-    MetaModelOpts.ExpDesign.UtilityFunction = 'DKL' #['DKL', 'DPP']
+    #MetaModelOpts.ExpDesign.UtilityFunction = 'DKL' #['DKL', 'DPP']
     
     # VarBasedOptDesign -> when data is not available
     # Only with Vornoi >>> 1)Entropy 2)EIGF, 3)ALM, 4)LOOCV
-    #MetaModelOpts.ExpDesign.UtilityFunction = 'LOOCV'
+    MetaModelOpts.ExpDesign.UtilityFunction = 'LOOCV'
     
     # alphabetic
     # 1)D-Opt (D-Optimality) 2)A-Opt (A-Optimality)