Initialize mutable Params after DAE discretization

pyomo/dae/misc.py

@@ -190,14 +190,6 @@ def update_contset_indexed_component(comp, expansion_map):
     if comp.ctype is Suffix:
         return
 
-    # Params indexed by a ContinuousSet should include an initialize
-    # and/or default rule which will be called automatically when the
-    # parameter value at a new point in the ContinuousSet is
-    # requested. Therefore, no special processing is required for
-    # Params.
-    if comp.ctype is Param:
-        return
-
     # Integral components are handled after every ContinuousSet has been
     # discretized. Import is deferred to here due to circular references.
     from pyomo.dae import Integral
@@ -228,6 +220,12 @@ def update_contset_indexed_component(comp, expansion_map):
                 # as Var components
                 expansion_map[comp] = _update_var
                 _update_var(comp)
+            elif comp.ctype is Param:
+                # Mutable Params with an initialize rule need to be updated
+                # after new ContinuousSet points are added. Immutable Params
+                # intentionally retain their existing behavior.
+                expansion_map[comp] = _update_param
+                _update_param(comp)
             elif comp.ctype == Constraint:
                 expansion_map[comp] = _update_constraint
                 _update_constraint(comp)
@@ -267,6 +265,40 @@ def _update_var(v):
         v.add(index)
 
 
+def _update_param(p):
+    """
+    Initialize new indices in a mutable Param after a ContinuousSet changes.
+
+    The Param initializer is only applied during component construction. A
+    DAE discretization can add points to a ContinuousSet after that, so the
+    normal construction path does not see the new indices. Defaults are
+    handled lazily by Param itself; this method only applies explicit
+    initialize rules to missing data.
+    """
+    if not p.mutable or p._rule is None or p._rule.contains_indices():
+        return
+    if not p.index_set().isfinite():
+        return
+
+    new_indices = [
+        index
+        for index in p.index_set()
+        if index not in p._data or p._data[index]._value is Param.NoValue
+    ]
+    if not new_indices:
+        return
+
+    block = p.parent_block()
+    rule = p._rule
+    if rule.constant():
+        value = rule(block, None)
+        for index in new_indices:
+            p._setitem_when_not_present(index, value)
+    else:
+        for index in new_indices:
+            p._setitem_when_not_present(index, rule(block, index))
+
+
 def _update_constraint(con):
     """
     This method will construct any additional indices in a constraint

pyomo/dae/tests/test_misc.py

@@ -387,6 +387,35 @@ def _con3(m, i, ti, ti2, j, k):
         self.assertTrue(value(m.con3[2, 0, 2, 1, 1].lower) is None)
         self.assertTrue(value(m.con3[3, 2, 3, 2, 2].upper) == 20)
 
+    def test_discretized_params_initialized_on_new_points(self):
+        for transformation, options in (
+            ('dae.finite_difference', dict(nfe=4)),
+            ('dae.collocation', dict(nfe=2, ncp=2)),
+        ):
+            m = ConcreteModel()
+            m.t = ContinuousSet(bounds=(0, 1))
+            m.s = Set(initialize=['a', 'b'])
+            m.p = Param(m.t, initialize=3, mutable=True)
+            m.p_indexed = Param(
+                m.t,
+                m.s,
+                initialize=lambda m, t, s: t + (1 if s == 'a' else 2),
+                mutable=True,
+            )
+            m.p_default = Param(m.t, m.s, default=5, mutable=True)
+            m.p_immutable = Param(m.t, initialize=7)
+
+            TransformationFactory(transformation).apply_to(m, **options)
+
+            new_point = next(t for t in m.t if t not in (0, 1))
+            self.assertEqual(m.p[new_point].value, 3)
+            self.assertAlmostEqual(m.p_indexed[new_point, 'a'].value, new_point + 1)
+            self.assertAlmostEqual(m.p_indexed[new_point, 'b'].value, new_point + 2)
+            self.assertEqual(m.p_default[new_point, 'a'].value, 5)
+            self.assertEqual(value(m.p_immutable[0]), 7)
+            with self.assertRaisesRegex(ValueError, 'undefined'):
+                value(m.p_immutable[new_point])
+
     # test update_contset_indexed_component method for Expression with
     # single index of the ContinuouSet
     def test_update_contset_indexed_component_expressions_single(self):