From 064ee411bf51622b7a076848c31541e4575c9238 Mon Sep 17 00:00:00 2001 From: mameikagou Date: Wed, 26 Aug 2026 06:56:32 +0800 Subject: [PATCH] Fix Param initialization after DAE discretization --- pyomo/dae/misc.py | 48 ++++++++++++++++++++++++++++++------ pyomo/dae/tests/test_misc.py | 29 ++++++++++++++++++++++ 2 files changed, 69 insertions(+), 8 deletions(-) diff --git a/pyomo/dae/misc.py b/pyomo/dae/misc.py index 0342dc81402..079150dc830 100644 --- a/pyomo/dae/misc.py +++ b/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 diff --git a/pyomo/dae/tests/test_misc.py b/pyomo/dae/tests/test_misc.py index c7a27e3a5ec..72055c04bff 100644 --- a/pyomo/dae/tests/test_misc.py +++ b/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):