diff --git a/sklearn/gaussian_process/_gpr.py b/sklearn/gaussian_process/_gpr.py
index 642ab4eef04c0..e5b3ca07045e4 100644
--- a/sklearn/gaussian_process/_gpr.py
+++ b/sklearn/gaussian_process/_gpr.py
@@ -196,10 +196,15 @@ def fit(self, X, y):
         if self.normalize_y:
             self._y_train_mean = np.mean(y, axis=0)
             self._y_train_std = np.std(y, axis=0)
-
+            # Assign a standard deviation of one to constant
+            # targets for avoiding division by zero errors
+            if np.iterable(self._y_train_std):
+                self._y_train_std[self._y_train_std == 0] = 1
+            else:
+                self._y_train_std = (
+                    self._y_train_std if self._y_train_std != 0 else 1)
             # Remove mean and make unit variance
             y = (y - self._y_train_mean) / self._y_train_std
-
         else:
             self._y_train_mean = np.zeros(1)
             self._y_train_std = 1
diff --git a/sklearn/gaussian_process/tests/test_gpr.py b/sklearn/gaussian_process/tests/test_gpr.py
index b3a4e7160cfdf..e580b4afc77fc 100644
--- a/sklearn/gaussian_process/tests/test_gpr.py
+++ b/sklearn/gaussian_process/tests/test_gpr.py
@@ -536,3 +536,17 @@ def test_bound_check_fixed_hyperparameter():
                         periodicity_bounds="fixed")  # seasonal component
     kernel = k1 + k2
     GaussianProcessRegressor(kernel=kernel).fit(X, y)
+
+
+@pytest.mark.parametrize('kernel', kernels)
+def test_constant_target(kernel):
+    # Test for constant target(s), the private _y_train_std attribute has
+    # a standard deviation of 1 instead of 0
+    gpr = GaussianProcessRegressor(kernel=kernel, normalize_y=True)
+    # create a constant targets of 2
+    y = np.full((1, X.shape[0]), 2).ravel()
+    expected_std = np.full((1, X.shape[0]), 1).ravel()
+    gpr.fit(X, y)
+    _y_train_std = gpr._y_train_std
+
+    assert_array_equal(_y_train_std, expected_std)