diff --git a/pyomo/repn/standard_repn.py b/pyomo/repn/standard_repn.py
index de83da61a01..b543f264e59 100644
--- a/pyomo/repn/standard_repn.py
+++ b/pyomo/repn/standard_repn.py
@@ -734,16 +734,29 @@ def _collect_pow(exp, multiplier, idMap, compute_values, verbose, quadratic):
             elif compute_values and len(res.linear) == 0:
                 return Results(constant=multiplier*res.constant**exponent)
             #
-            # If there is one linear term, then we compute the quadratic expression for it.
+            # If the base is linear, then we compute the quadratic expression for it.
             #
-            elif len(res.linear) == 1:
-                key, coef = res.linear.popitem()
-                var = idMap[key]
+            else:
                 ans = Results()
-                if not (res.constant.__class__ in native_numeric_types and res.constant == 0):
+                has_constant = (res.constant.__class__
+                                not in native_numeric_types
+                                or res.constant != 0)
+                if has_constant:
                     ans.constant = multiplier*res.constant*res.constant
-                    ans.linear[key] = 2*multiplier*coef*res.constant
-                ans.quadratic[key,key] = multiplier*coef*coef
+
+                # this is reversed since we want to pop off the end for efficiency
+                # and the quadratic terms have a convention that the indexing tuple
+                # of key1, key2 is such that key1 <= key2
+                keys = sorted(res.linear.keys(), reverse=True)
+                while len(keys) > 0:
+                    key1 = keys.pop()
+                    coef1 = res.linear[key1]
+                    if has_constant:
+                        ans.linear[key1] = 2*multiplier*coef1*res.constant
+                    ans.quadratic[key1,key1] = multiplier*coef1*coef1
+                    for key2 in keys:
+                        coef2 = res.linear[key2]
+                        ans.quadratic[key1,key2] = 2*multiplier*coef1*coef2
                 return ans
 
     #
diff --git a/pyomo/repn/tests/test_standard.py b/pyomo/repn/tests/test_standard.py
index f9a8a60707f..fd12f26de7a 100644
--- a/pyomo/repn/tests/test_standard.py
+++ b/pyomo/repn/tests/test_standard.py
@@ -3141,19 +3141,25 @@ def test_pow5(self):
         rep = generate_standard_repn(e, compute_values=False, quadratic=True)
         #
         self.assertFalse( rep.is_fixed() )
-        self.assertEqual( rep.polynomial_degree(), None )
+        self.assertEqual( rep.polynomial_degree(), 2 )
         self.assertFalse( rep.is_constant() )
         self.assertFalse( rep.is_linear() )
-        self.assertFalse( rep.is_quadratic() )
+        self.assertTrue( rep.is_quadratic() )
         self.assertTrue( rep.is_nonlinear() )
         #
         self.assertTrue(len(rep.linear_vars) == 0)
         self.assertTrue(len(rep.linear_coefs) == 0)
-        self.assertTrue(len(rep.quadratic_vars) == 0)
-        self.assertTrue(len(rep.quadratic_coefs) == 0)
-        self.assertFalse(rep.nonlinear_expr is None)
-        self.assertTrue(len(rep.nonlinear_vars) == 2)
-        baseline = { }
+        self.assertTrue(len(rep.quadratic_vars) == 3)
+        self.assertTrue(len(rep.quadratic_coefs) == 3)
+        self.assertTrue(rep.nonlinear_expr is None)
+        self.assertTrue(len(rep.nonlinear_vars) == 0)
+        baseline = { (id(m.a), id(m.a)): 1,
+                     (id(m.b), id(m.b)): 1 }
+        if id(m.a) < id(m.b):
+            baseline[id(m.a), id(m.b)] = 2
+        else:
+            baseline[id(m.b), id(m.a)] = 2
+
         self.assertEqual(baseline, repn_to_dict(rep))
 
         e = (m.a+3)**2
@@ -3258,6 +3264,51 @@ def test_pow6(self):
         baseline = { None:8 }
         self.assertEqual(baseline, repn_to_dict(rep))
 
+    def test_pow_of_lin_sum(self):
+        m = ConcreteModel()
+        m.x = Var(range(4))
+        e = sum(x for x in m.x.values())**2
+
+        rep = generate_standard_repn(e, compute_values=False, quadratic=False)
+        #
+        self.assertFalse( rep.is_fixed() )
+        self.assertEqual( rep.polynomial_degree(), None )
+        self.assertFalse( rep.is_constant() )
+        self.assertFalse( rep.is_linear() )
+        self.assertFalse( rep.is_quadratic() )
+        self.assertTrue( rep.is_nonlinear() )
+        #
+        self.assertTrue(len(rep.linear_vars) == 0)
+        self.assertTrue(len(rep.linear_coefs) == 0)
+        self.assertTrue(len(rep.quadratic_vars) == 0)
+        self.assertTrue(len(rep.quadratic_coefs) == 0)
+        self.assertFalse(rep.nonlinear_expr is None)
+        self.assertTrue(len(rep.nonlinear_vars) == 4)
+        baseline = { }
+        self.assertEqual(baseline, repn_to_dict(rep))
+
+        rep = generate_standard_repn(e, compute_values=False, quadratic=True)
+        #
+        self.assertFalse( rep.is_fixed() )
+        self.assertEqual( rep.polynomial_degree(), 2 )
+        self.assertFalse( rep.is_constant() )
+        self.assertFalse( rep.is_linear() )
+        self.assertTrue( rep.is_quadratic() )
+        self.assertTrue( rep.is_nonlinear() )
+        #
+        self.assertTrue(len(rep.linear_vars) == 0)
+        self.assertTrue(len(rep.linear_coefs) == 0)
+        self.assertTrue(len(rep.quadratic_vars) == 10)
+        self.assertTrue(len(rep.quadratic_coefs) == 10)
+        self.assertTrue(rep.nonlinear_expr is None)
+        self.assertTrue(len(rep.nonlinear_vars) == 0)
+        baseline = {(id(i), id(j)): 2
+                    for i in m.x.values()
+                    for j in m.x.values()
+                    if id(i) < id(j)}
+        baseline.update({(id(i), id(i)): 1 for i in m.x.values()})
+        self.assertEqual(baseline, repn_to_dict(rep))
+
     def test_fixed_exponent(self):
         m = ConcreteModel()
         m.x = Var()