Port the Pyomo extension to Pyomo 2 / Symbolics 7 (#4888)

* Port the Pyomo extension to Pyomo 2 / Symbolics 7

MTKPyomoDynamicOptExt could not load at all: its compat was pinned to
Pyomo 0.1.0, and every registered Pyomo version up to 1.0.1 declared
`Symbolics = "6.39.1 - 6"` while MTKBase requires Symbolics 7.32.1. The
extension also called `Symbolics.symstruct` and `Symbolics.Struct` directly,
both removed in Symbolics 7.

* Compat bumped to Pyomo 2.0.1.
* `Symbolics.symstruct(ConcreteModel)` -> `SymbolicConcreteModel`. Symbolics 7
  uses the Julia type itself as the symtype.
* `pysym_getproperty` / `_getproperty` deleted here; they now live in Pyomo.jl,
  which can test them against a real model.
* Symbolic indexing goes through `Pyomo.pyomo_getindex`. Symbolics 7 moved the
  symtype out of `BasicSymbolic`'s type parameter, so the old
  `BasicSymbolic{Struct{PyomoVar}}` `getindex` dispatch is not expressible; the
  `Symbolics.value` calls that existed only to reach it are gone too.
* `Symbolics.substitute(..., fold = false)` -> `fold = Val(false)`. SymbolicUtils
  takes `fold` as a `Val`; the `Bool` form throws on every `add_constraint!` and
  `set_objective!`. `Val` is valid at the declared SymbolicUtils floor (4.35.3).

`using Pyomo` is uncommented in the test file and Pyomo added to the
optimization test environment, where it was not a dependency. Every Pyomo
assertion sits behind `if @isdefined(Pyomo)`, so none of them have run.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016GYQeAiC2x9nJ7Mmwwxhe2

* Require Pyomo 2.0.2

2.0.1 shipped without the promote_shape, isinteger and operator fixes; CI on
this branch failed against it with

    ArgumentError: Invalid shapes for cos: (SymbolicUtils.Unknown(-1),).

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016GYQeAiC2x9nJ7Mmwwxhe2

---------

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: 3 people

lib/ModelingToolkitBase/Project.toml

@@ -164,7 +164,7 @@ Pkg = "1"
 PreallocationTools = "0.4.27, 1"
 PrecompileTools = "1.2.1"
 Printf = "1"
-Pyomo = "0.1.0"
+Pyomo = "2.0.2"
 REPL = "1"
 Random = "1"
 ReadOnlyDicts = "1.0.0"

lib/ModelingToolkitBase/ext/MTKPyomoDynamicOptExt.jl

@@ -65,7 +65,7 @@ struct PyomoDynamicOptModel{T}
     dummy_sym::Union{Num, Symbolics.BasicSymbolic}
 
     function PyomoDynamicOptModel(model, U, V, P, tₛ, is_free_final, tsteps)
-        @variables MODEL_SYM::Symbolics.symstruct(ConcreteModel) T_SYM DUMMY_SYM
+        @variables MODEL_SYM::SymbolicConcreteModel T_SYM DUMMY_SYM
         model.dU = dae.DerivativeVar(U, wrt = model.t, initialize = 0)
         return new{typeof(P)}(
             model, U, V, P, tₛ, is_free_final, tsteps, nothing,
@@ -94,15 +94,6 @@ struct PyomoDynamicOptProblem{uType, tType, isinplace, P, F, K} <:
     end
 end
 
-function pysym_getproperty(s::Union{Num, SymbolicT}, name::Symbol)
-    return Symbolics.wrap(
-        SymbolicUtils.term(
-            _getproperty, Symbolics.unwrap(s), Val{name}(), type = Symbolics.Struct{PyomoVar}
-        )
-    )
-end
-_getproperty(s, name::Val{fieldname}) where {fieldname} = getproperty(s, fieldname)
-
 function MTK.PyomoDynamicOptProblem(
         sys::System, op, tspan;
         dt = nothing, steps = nothing, tune_parameters = false,
@@ -165,7 +156,7 @@ function MTK.add_constraint!(pmodel::PyomoDynamicOptModel, cons; n_idxs = 1)
         cons.lhs - cons.rhs ≤ 0
     end
     expr = Symbolics.substitute(
-        Symbolics.unwrap(expr), SPECIAL_FUNCTIONS_DICT, fold = false
+        Symbolics.unwrap(expr), SPECIAL_FUNCTIONS_DICT, fold = Val(false)
     )
 
     cons_sym = Symbol("cons", hash(cons))
@@ -192,8 +183,8 @@ function MTK.set_variable_bounds!(m::PyomoDynamicOptModel, sys, pmap, tf, tunabl
         MTK.add_constraint!(m, var ≲ hi)
     end
     for (i, (lo, hi)) in param_bounds
-        P_sym = Symbolics.value(pysym_getproperty(m.model_sym, :P))
-        p_var = P_sym[i]
+        P_sym = pysym_getproperty(m.model_sym, :P)
+        p_var = pyomo_getindex(P_sym, i)
         MTK.add_constraint!(m, p_var ≳ lo)
         MTK.add_constraint!(m, p_var ≲ hi)
     end
@@ -206,7 +197,7 @@ end
 
 function MTK.set_objective!(pmodel::PyomoDynamicOptModel, expr)
     @unpack model, model_sym, t_sym, dummy_sym = pmodel
-    expr = Symbolics.substitute(expr, SPECIAL_FUNCTIONS_DICT, fold = false)
+    expr = Symbolics.substitute(expr, SPECIAL_FUNCTIONS_DICT, fold = Val(false))
     return if SU.query(isequal(Symbolics.unwrap(t_sym)), expr)
         f = eval(Symbolics.build_function(expr, model_sym, t_sym))
         model.obj = pyomo.Objective(model.t, rule = Pyomo.pyfunc(f))
@@ -243,13 +234,13 @@ function MTK.lowered_integral(m::PyomoDynamicOptModel, arg, lo, hi)
 end
 
 function MTK.lowered_derivative(m::PyomoDynamicOptModel, i)
-    mdU = Symbolics.value(pysym_getproperty(m.model_sym, :dU))
-    return Symbolics.unwrap(mdU[i, m.t_sym])
+    mdU = pysym_getproperty(m.model_sym, :dU)
+    return Symbolics.unwrap(pyomo_getindex(mdU, i, m.t_sym))
 end
 
 function MTK.lowered_var(m::PyomoDynamicOptModel, uv, i, t)
-    X = Symbolics.value(pysym_getproperty(m.model_sym, uv))
-    var = t isa Union{Num, SymbolicT} ? X[i, m.t_sym] : X[i, t]
+    X = pysym_getproperty(m.model_sym, uv)
+    var = t isa Union{Num, SymbolicT} ? pyomo_getindex(X, i, m.t_sym) : pyomo_getindex(X, i, t)
     return Symbolics.unwrap(var)
 end
 
@@ -258,9 +249,9 @@ end
 function MTK.get_param_for_pmap(m::ConcreteModel, P::PyomoVar, i)
     # Create a symbolic variable that will be used in the pmap
     # The actual PyomoVar will be accessed via the symbolic representation
-    @variables MODEL_SYM::Symbolics.symstruct(ConcreteModel)
-    P_sym = Symbolics.value(pysym_getproperty(MODEL_SYM, :P))
-    return Symbolics.unwrap(P_sym[i])
+    @variables MODEL_SYM::SymbolicConcreteModel
+    P_sym = pysym_getproperty(MODEL_SYM, :P)
+    return Symbolics.unwrap(pyomo_getindex(P_sym, i))
 end
 
 MTK.needs_individual_tunables(m::ConcreteModel) = true

lib/ModelingToolkitBase/test/optimization/Project.toml

@@ -19,6 +19,7 @@ OrdinaryDiffEqImplicitTableaus = "75f66a49-58fc-43e3-9173-2340726368f7"
 OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+Pyomo = "0e8e1daf-01b5-4eba-a626-3897743a3816"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
@@ -35,5 +36,6 @@ CasADi = "1.0.7"
 DataInterpolations = "8.8, 9"
 OrdinaryDiffEqExplicitTableaus = "2"
 OrdinaryDiffEqImplicitTableaus = "2"
+Pyomo = "2.0.2"
 SafeTestsets = "0.1, 1"
 SciMLTesting = "1, 2.1"

lib/ModelingToolkitBase/test/optimization/dynamic_optimization.jl

@@ -10,7 +10,7 @@ using Ipopt
 using DataInterpolations
 using DataInterpolations: ConstantInterpolation
 using CasADi
-# using Pyomo
+using Pyomo
 using Test
 
 import DiffEqBase: solve