diff --git a/deerlab/fitparamodel.py b/deerlab/fitparamodel.py
index e54b6a7b8..2e9266e8d 100644
--- a/deerlab/fitparamodel.py
+++ b/deerlab/fitparamodel.py
@@ -202,7 +202,7 @@ def lsqresiduals(p):
         sol = least_squares(lsqresiduals ,par0, bounds=(lb,ub), max_nfev=int(maxiter), ftol=tol, method='dogbox')
         sols.append(sol)
         parfits.append(sol.x)
-        fvals.append(sol.cost)        
+        fvals.append(2*sol.cost) # least_squares uses 0.5*sum(residual**2)          
 
     # Find global minimum from multiple runs
     globmin = np.argmin(fvals)
@@ -255,7 +255,7 @@ def lsqresiduals(p):
     if Nsignals==1: 
         stats = stats[0]
         scales = scales[0]
-
+        fvals = fvals[0]
     # Get plot function
     plotfcn = lambda: _plot(Vsubsets,V,Vfit)
 
diff --git a/deerlab/snlls.py b/deerlab/snlls.py
index 06c8798f4..9cdcb3ada 100644
--- a/deerlab/snlls.py
+++ b/deerlab/snlls.py
@@ -353,7 +353,7 @@ def ResidualsFcn(p):
         sol = least_squares(ResidualsFcn, par0, bounds=(lb, ub), max_nfev=int(nonlin_maxiter), ftol=nonlin_tol)
         nonlinfits.append(sol.x)
         linfits.append(linfit)
-        fvals.append(sol.cost)
+        fvals.append(2*sol.cost) # least_squares uses 0.5*sum(residual**2)          
         sols.append(sol)
     # Find global minimum from multiple runs
     globmin = np.argmin(fvals)
@@ -418,6 +418,7 @@ def ci(coverage,ptype='full'):
         stats.append(goodness_of_fit(y[subset], yfit[subset], Ndof))
     if len(stats) == 1: 
         stats = stats[0]
+        fvals = fvals[0]
 
     # Display function
     plotfcn = lambda: _plot(subsets,y,yfit)
diff --git a/test/test_fitmultimodel.py b/test/test_fitmultimodel.py
index fbfc43dea..d8c7198a3 100644
--- a/test/test_fitmultimodel.py
+++ b/test/test_fitmultimodel.py
@@ -254,3 +254,18 @@ def test_globalfit_scales():
     assert max(abs(np.asarray(scales)/np.asarray(fit.scale) - 1)) < 1e-2 
 #=======================================================================
 
+def test_cost_value():
+#=======================================================================
+    "Check that the fit of a multi-Gauss model works"
+        
+    r = np.linspace(2,6,300)
+    t = np.linspace(-0.5,6,500)
+    K = dipolarkernel(t,r)
+    parin = [4, 0.05, 0.4, 4, 0.4, 0.4, 3, 0.15, 0.2]
+    P = dd_gauss3(r,parin)
+    V = K@P
+
+    fit = fitmultimodel(V,K,r,dd_gauss,3,'aicc', uqanalysis=False)
+    
+    assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
+#=======================================================================
diff --git a/test/test_fitparamodel.py b/test/test_fitparamodel.py
index 74b21e851..c6e5ffb7e 100644
--- a/test/test_fitparamodel.py
+++ b/test/test_fitparamodel.py
@@ -224,4 +224,26 @@ def Vmodel(par):
     fit = fitparamodel([V1,V2],Vmodel,par0,lb,ub,weights=[1,1])
 
     assert max(abs(np.asarray(scales)/np.asarray(fit.scale) - 1)) < 1e-2 
+#============================================================
+
+
+def test_cost_value():
+#============================================================
+    "Check that the cost value is properly returned"
+
+    t = np.linspace(0,3,300)
+    r = np.linspace(2,5,200)
+    par = np.array([3,0.2])
+    P = dd_gauss(r,par)
+
+    K = dipolarkernel(t,r)
+    V = K@P
+
+    par0 = [5, 0.5]
+    lb = [1, 0.1]
+    ub = [20, 1]
+    model = lambda p: K@dd_gauss(r,p)
+    fit = fitparamodel(V,model,par0,lb,ub)
+
+    assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
 #============================================================
\ No newline at end of file
diff --git a/test/test_fitregmodel.py b/test/test_fitregmodel.py
index 43e2356f7..810e59f5b 100644
--- a/test/test_fitregmodel.py
+++ b/test/test_fitregmodel.py
@@ -390,3 +390,18 @@ def test_globalfit_scales():
 
     assert max(abs(np.asarray(scales)/np.asarray(fit.scale) - 1)) < 1e-2 
 #============================================================
+
+def test_cost_value():
+#============================================================
+    "Check that the cost value is properly returned"
+
+    np.random.seed(1)
+    t = np.linspace(-2,4,300)
+    r = np.linspace(2,6,100)
+    P = dd_gauss(r,[3,0.2])
+    K = dipolarkernel(t,r)
+    V = K@P + whitegaussnoise(t,0.01)
+    fit = fitregmodel(V,K,r)
+
+    assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
+#============================================================
diff --git a/test/test_fitsignal.py b/test/test_fitsignal.py
index 1b50476a7..ccf1c7d78 100644
--- a/test/test_fitsignal.py
+++ b/test/test_fitsignal.py
@@ -508,4 +508,22 @@ def test_V_scale_regularized():
     fit = fitsignal(V,t,r,'P',None,None,uqanalysis=False)
 
     assert max(abs(1 - V/fit.V)) < 1e-4
-# ======================================================================
\ No newline at end of file
+# ======================================================================
+
+def test_cost_value():
+# ======================================================================
+    "Check that starts values can be correctly specified"
+
+    t = np.linspace(0,5,100)
+    r = np.linspace(2,6,150)
+    P = dd_gauss(r,[4.5, 0.25])
+
+    Bmodel = lambda t: bg_exp(t,0.4)
+    K = dipolarkernel(t,r,0.4,Bmodel)
+    V = K@P
+
+    fit = fitsignal(V,t,r,'P',bg_exp,ex_4pdeer,uqanalysis=False)
+
+    assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
+# ======================================================================
+
diff --git a/test/test_snlls.py b/test/test_snlls.py
index b708c6897..cfa44b766 100644
--- a/test/test_snlls.py
+++ b/test/test_snlls.py
@@ -352,4 +352,33 @@ def test_reg_huber():
     "Check SNNLS using Tikhonov regularization works"
 
     assert_reg_type(regtype='huber')
-#=======================================================================
\ No newline at end of file
+#=======================================================================
+
+
+def test_cost_value():
+#============================================================
+    "Check that the cost value is properly returned"
+
+    # Prepare test data
+    r = np.linspace(1,8,80)
+    t = np.linspace(0,4,200)
+    lam = 0.25
+    K = dipolarkernel(t,r,lam)
+    parin = [3.5, 0.4, 0.6, 4.5, 0.5, 0.4]
+    P = dd_gauss2(r,parin)
+    V = K@P
+
+    # Non-linear parameters
+    # nlpar = [lam]
+    nlpar0 = 0.2
+    lb = 0
+    ub = 1
+    # Linear parameters: non-negativity
+    lbl = np.zeros(len(r))
+    ubl = np.full(len(r), np.inf)
+    # Separable LSQ fit
+
+    fit = snlls(V,lambda lam: dipolarkernel(t,r,lam),nlpar0,lb,ub,lbl,ubl)
+
+    assert isinstance(fit.cost,float) and np.round(fit.cost/np.sum(fit.residuals**2),5)==1
+#============================================================
\ No newline at end of file