In the Method of Morris, result metric sigma is not computed correctly
The result metric sigma should be an unbiased estimator of the sample variation.  The current way of computing the standard deviation (`np.std(ee)`) is incorrect, and should be changed to `np.std(ee, ddof=1)` 

diff --git a/SALib/analyze/morris.py b/SALib/analyze/morris.py
index 70d05c5..c83eace 100755
--- a/SALib/analyze/morris.py
+++ b/SALib/analyze/morris.py
@@ -46,7 +46,7 @@ def analyze(problem, X, Y,
               for k in ['names', 'mu', 'mu_star', 'sigma', 'mu_star_conf'])
     Si['mu'] = np.average(ee, 1)
     Si['mu_star'] = np.average(np.abs(ee), 1)
-    Si['sigma'] = np.std(ee, 1)
+    Si['sigma'] = np.std(ee, axis=1, ddof=1)
     Si['names'] = problem['names']
 
     for j in range(num_vars):
@@ -61,20 +61,23 @@ def analyze(problem, X, Y,
                     'sigma'][j], Si['mu_star'][j], Si['mu_star_conf'][j]))
         return Si
     elif groups is not None:
+        Si['sigma'] = np.std(np.abs(ee), axis=1, ddof=1)
         # if there are groups, then the elementary effects returned need to be
         # computed over the groups of variables, rather than the individual variables
         Si_grouped = dict((k, [None] * num_vars)
                 for k in ['mu_star', 'mu_star_conf'])
-        Si_grouped['mu_star'] = compute_grouped_mu_star(Si['mu_star'], groups)
-        Si_grouped['mu_star_conf'] = compute_grouped_mu_star(Si['mu_star_conf'],
+        Si_grouped['mu_star'] = compute_grouped_metric(Si['mu_star'], groups)
+        Si_grouped['sigma'] = compute_grouped_sigma(Si, groups)
+        Si_grouped['mu_star_conf'] = compute_grouped_metric(Si['mu_star_conf'],
                                                              groups)
         Si_grouped['names'] = unique_group_names
 
         if print_to_console:
-            print("Parameter Mu_Star Mu_Star_Conf")
+            print("Parameter Mu_Star Sigma Mu_Star_Conf")
             for j in list(range(number_of_groups)):
-                print("%s %f %f" % (Si_grouped['names'][j],
+                print("%s %f %f %f" % (Si_grouped['names'][j],
                                     Si_grouped['mu_star'][j],
+                                    Si_grouped['sigma'][j],
                                     Si_grouped['mu_star_conf'][j]))
 
         return Si_grouped
@@ -82,12 +85,32 @@ def analyze(problem, X, Y,
         raise RuntimeError("Could not determine which parameters should be returned")
 
 
-def compute_grouped_mu_star(mu_star_ungrouped, group_matrix):
+def compute_grouped_metric(ungrouped_metric, group_matrix):
 
     group_matrix = np.array(group_matrix)
-    mu_star_grouped = np.divide(np.dot(mu_star_ungrouped, group_matrix), np.sum(group_matrix, 0))
+    grouped_metric = np.divide(np.dot(ungrouped_metric, group_matrix), np.sum(group_matrix, 0))
 
-    return mu_star_grouped.T
+    return grouped_metric.T
+
+
+def compute_grouped_sigma(Si, group_matrix):
+    
+    group_matrix = np.array(group_matrix, dtype=np.bool)
+
+    mu_star_masked = np.ma.masked_array(Si['mu_star'] * group_matrix.T, mask=(group_matrix^1).T)
+    sigma_masked = np.ma.masked_array(Si['sigma'] * group_matrix.T, mask=(group_matrix^1).T)**2
+    
+    mean_of_variances_of_ee = np.ma.mean(sigma_masked, axis=1)
+#     variance_of_means = np.var(Si['mu_star'] * group_matrix.T, axis=1)
+# ``var = mean(abs(x - x.mean())**2)``
+
+    variance_of_means_of_ee = np.ma.var(mu_star_masked, axis=1, ddof=1)
+
+    grouped_sigma = np.sqrt(mean_of_variances_of_ee + variance_of_means_of_ee)
+
+    grouped_sigma = np.where(np.sum(group_matrix, axis=0)==1, np.sqrt(np.sum(sigma_masked, axis=1)), grouped_sigma)
+
+    return grouped_sigma.T
 
 
 def get_increased_values(op_vec, up, lo):
@@ -154,7 +177,7 @@ def compute_mu_star_confidence(ee, num_trajectories, num_resamples, conf_level):
     if conf_level < 0 or conf_level > 1:
         raise ValueError("Confidence level must be between 0-1.")
 
-    resample_index = np.random.randint(len(ee), size=(num_resamples, num_trajectories)) 
+    resample_index = np.random.randint(len(ee), size=(num_resamples, num_trajectories))
     ee_resampled = ee[resample_index]
     # Compute average of the absolute values over each of the resamples
     mu_star_resampled = np.average(np.abs(ee_resampled), axis=1)
diff --git a/SALib/sample/optimal_trajectories.py b/SALib/sample/optimal_trajectories.py
index 0d04c33..1ceebb2 100644
--- a/SALib/sample/optimal_trajectories.py
+++ b/SALib/sample/optimal_trajectories.py
@@ -34,7 +34,6 @@ def model(N, k_choices, distance_matrix):
 
     m = Model("distance1")
     I = range(N)
-    big_M = k_choices + 1
 
     distance_matrix = distance_matrix / distance_matrix.max()
 
@@ -49,23 +48,20 @@ def model(N, k_choices, distance_matrix):
 
     m.setObjective(quicksum([x[i, j] * dm[j][i] for i in I for j in range(i + 1, N)]))
 
-    m.addConstr(quicksum([x[i, j] for i in I for j in range(i + 1, N)]) == nchoosek(k_choices, 2), "All")
-
-    # Finally, each combination may only appear three times in the combination list
-    for i in I:
-        m.addConstr(quicksum(x[i, j] for j in range(i + 1, N)) + quicksum(x[k, i] for k in range(0, i)) - (y[i] * big_M),
-                    '<=',
-                    (k_choices - 1),
-                    "a:Only k-1 scores in any row/column for %s" % i)
-        m.addConstr(quicksum(x[i, j] for j in range(i + 1, N)) + quicksum(x[k, i] for k in range(0, i)) + (y[i] * big_M),
-                    '>=',
-                    (k_choices - 1),
-                    "b:Only k-1 scores in any row/column for %s" % i)
-
-    m.addConstr(quicksum(y[i] for i in I), "==", N - k_choices, name="Only %s hold" % (N - k_choices))
+    # Add constraints to the model     
+    m.addConstr(quicksum([y[i] for i in I]) <= k_choices, "27")
+    
+    for i in I: 
+        for j in range(i + 1, N):
+            m.addConstr(x[i, j] <= y[i], "28-%s-%s" % (i, j))
+            m.addConstr(x[i, j] <= y[j], "29-%s-%s" % (i, j))
+            m.addConstr(y[i] + y[j] <= 1 + x[i, j], "30-%s-%s" % (i, j))
+    
+    m.addConstr(quicksum([x[i, j] for i in I for j in range(i + 1, N)]) <= nchoosek(k_choices, 2), "Cut_1")
     m.update()
     return m
 
+
 @requires_gurobipy(_has_gurobi)
 def return_max_combo(input_data, N, num_params, k_choices, groups=None):
     
@@ -80,7 +76,7 @@ def return_max_combo(input_data, N, num_params, k_choices, groups=None):
     m.params.Threads = 0
     # m.params.MIPGap=0.03
 
-    # m.write("model.lp")
+    m.write("model.lp")
     m.ModelSense = GRB.MAXIMIZE
     m.optimize()
     if m.Status == GRB.OPTIMAL:
diff --git a/SALib/tests/test_analyze_morris.py b/SALib/tests/test_analyze_morris.py
index d7a6517..f49f555 100644
--- a/SALib/tests/test_analyze_morris.py
+++ b/SALib/tests/test_analyze_morris.py
@@ -6,12 +6,13 @@ from numpy.testing import assert_allclose, assert_equal
 
 import numpy as np
 
-from ..analyze.morris import analyze, \\
+from SALib.analyze.morris import analyze, \\
     compute_mu_star_confidence, \\
     compute_elementary_effects, \\
     get_increased_values, \\
     get_decreased_values, \\
-    compute_grouped_mu_star
+    compute_grouped_metric, \\
+    compute_grouped_sigma
 
 
 def test_compute_mu_star_confidence():
@@ -32,6 +33,8 @@ def test_compute_mu_star_confidence():
 def test_analysis_of_morris_results():
     '''
     Tests a one-dimensional vector of results
+    
+    Taken from the solution to Exercise 4 (p.138) in Saltelli (2008).
     '''
     model_input = np.array([[0, 1. / 3], [0, 1], [2. / 3, 1],
                              [0, 1. / 3], [2. / 3, 1. / 3], [2. / 3, 1],
@@ -41,7 +44,9 @@ def test_analysis_of_morris_results():
                              [1. / 3, 2. / 3], [1. / 3, 0], [1, 0]],
                             dtype=np.float)
 
-    model_output = np.array([0.97, 0.71, 2.39, 0.97, 2.3, 2.39, 1.87, 2.40, 0.87, 2.15, 1.71, 1.54, 2.15, 2.17, 1.54, 2.2, 1.87, 1.0],
+    model_output = np.array([0.97, 0.71, 2.39, 0.97, 2.30, 2.39, 
+                             1.87, 2.40, 0.87, 2.15, 1.71, 1.54, 
+                             2.15, 2.17, 1.54, 2.20, 1.87, 1.0],
                        dtype=np.float)
 
     problem = {
@@ -55,16 +60,19 @@ def test_analysis_of_morris_results():
                  num_resamples=1000,
                  conf_level=0.95,
                  print_to_console=False)
-    ee = np.array([[2.52, 2.01, 2.30, -0.66, -0.93, -1.30],
-                   [-0.39, 0.13, 0.80, 0.25, -0.02, 0.51]])
-    desired_mu = np.average(ee, 1)
-    assert_allclose(Si['mu'], desired_mu, rtol=1e-1)
-    desired_mu_star = np.average(np.abs(ee), 1)
-    assert_allclose(Si['mu_star'], desired_mu_star, rtol=1e-2)
-    desired_sigma = np.std(ee, 1)
-    assert_allclose(Si['sigma'], desired_sigma, rtol=1e-2)
+
+    desired_mu = np.array([0.66, 0.21])
+    assert_allclose(Si['mu'], desired_mu, rtol=1e-1, 
+                    err_msg="The values for mu are incorrect")
+    desired_mu_star = np.array([1.62, 0.35])
+    assert_allclose(Si['mu_star'], desired_mu_star, rtol=1e-2, 
+                    err_msg="The values for mu star are incorrect")
+    desired_sigma = np.array([1.79, 0.41])
+    assert_allclose(Si['sigma'], desired_sigma, rtol=1e-2, 
+                    err_msg="The values for sigma are incorrect")
     desired_names = ['Test 1', 'Test 2']
-    assert_equal(Si['names'], desired_names)
+    assert_equal(Si['names'], desired_names, 
+                 err_msg="The values for names are incorrect")
 
 
 @raises(ValueError)
@@ -216,6 +224,28 @@ def test_compute_grouped_mu_star():
                    [-2.00, 0.13, -0.80, 0.25, -0.02, 0.51],
                    [2.00, -0.13, 0.80, -0.25, 0.02, -0.51]])
     mu_star = np.average(np.abs(ee), 1)
-    actual = compute_grouped_mu_star(mu_star, group_matrix)
+    actual = compute_grouped_metric(mu_star, group_matrix)
     desired = np.array([1.62, 0.62], dtype=np.float)
     assert_allclose(actual, desired, rtol=1e-1)
+    
+    
+def test_compute_grouped_sigma():
+    '''
+    Computes sigma for 3 variables grouped into 2 groups
+    There are six trajectories.
+    '''
+    group_matrix = np.matrix('1,0;0,1;0,1', dtype=np.int)
+    ee = np.array([[2.52, 2.01, 2.30, -0.66, -0.93, -1.30],
+                   [-2.00, 0.13, -0.80, 0.25, -0.02, 0.51],
+                   [2.00, -0.13, 0.80, -0.25, 0.02, -0.51]])
+    mu_star = np.average(np.abs(ee), 1)
+    Si = {}
+    Si['mu_star'] = mu_star
+    Si['sigma'] = np.std(np.abs(ee), ddof=1, axis=1)
+    
+    actual = compute_grouped_sigma(Si, group_matrix)
+    desired = np.array([0.7648, 0.7329], dtype=np.float64)
+    assert_allclose(actual, desired, rtol=1e-4)
+    
+    
+
['SALib/sample/optimal_trajectories.py', 'SALib/tests/test_analyze_morris.py', 'SALib/analyze/morris.py']