Pyomo Simulator does not work with Expression

[damdaepark]

Summary

Steps to reproduce the issue

from pyomo.environ import *
from pyomo.dae import *
from pyomo.dae.simulator import Simulator

def create_model():
    m = ConcreteModel()

    m.t = ContinuousSet(bounds=(0.0, 10.0))

    m.b = Param(initialize=0.25)
    m.c = Param(initialize=5.0)

    m.omega = Var(m.t)
    m.theta = Var(m.t)

    m.domegadt = DerivativeVar(m.omega, wrt=m.t)
    m.dthetadt = DerivativeVar(m.theta, wrt=m.t)

    # Setting the initial conditions
    m.omega[0] = 0.0
    m.theta[0] = 3.14 - 0.1

    m.eq1 = Expression(m.t, rule=lambda m, t: -m.b * m.omega[t] - m.c * sin(m.theta[t]))

    def _diffeq1(m, t):
        return m.domegadt[t] == m.eq1[t]
    m.diffeq1 = Constraint(m.t, rule=_diffeq1)

    def _diffeq2(m, t):
        return m.dthetadt[t] == m.omega[t]
    m.diffeq2 = Constraint(m.t, rule=_diffeq2)

    return m

def simulate_model(m):
    if False:
        # Simulate the model using casadi
        sim = Simulator(m, package='casadi')
        tsim, profiles = sim.simulate(numpoints=100, integrator='cvodes')
    else:
        # Simulate the model using scipy
        sim = Simulator(m, package='scipy')
        tsim, profiles = sim.simulate(numpoints=100, integrator='vode')

    # Discretize model using Orthogonal Collocation
    discretizer = TransformationFactory('dae.collocation')
    discretizer.apply_to(m, nfe=8, ncp=5)

    # Initialize the discretized model using the simulator profiles
    sim.initialize_model()

    return sim, tsim, profiles

def plot_result(m, sim, tsim, profiles):
    import matplotlib.pyplot as plt

    time = list(m.t)
    omega = [value(m.omega[t]) for t in m.t]
    theta = [value(m.theta[t]) for t in m.t]

    varorder = sim.get_variable_order()

    for idx, v in enumerate(varorder):
        plt.plot(tsim, profiles[:, idx], label=v)
    plt.plot(time, omega, 'o', label='omega interp')
    plt.plot(time, theta, 'o', label='theta interp')
    plt.xlabel('t')
    plt.legend(loc='best')
    plt.show()

model = create_model()
sim, tsim, profiles = simulate_model(model)
plot_result(model, sim, tsim, profiles)...

Error Message

capi_return is NULL
Call-back cb_f_in_dvode__user__routines failed.
Traceback (most recent call last):
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\runpy.py", line 194, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "c:\Users\Damdae\.vscode\extensions\ms-python.python-2022.0.1786462952\pythonFiles\lib\python\debugpy\__main__.py", line 45, in <module>
    cli.main()
  File "c:\Users\Damdae\.vscode\extensions\ms-python.python-2022.0.1786462952\pythonFiles\lib\python\debugpy/..\debugpy\server\cli.py", line 444, in main    
    run()
  File "c:\Users\Damdae\.vscode\extensions\ms-python.python-2022.0.1786462952\pythonFiles\lib\python\debugpy/..\debugpy\server\cli.py", line 285, in run_file
    runpy.run_path(target_as_str, run_name=compat.force_str("__main__"))
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\runpy.py", line 265, in run_path
    return _run_module_code(code, init_globals, run_name,
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\runpy.py", line 97, in _run_module_code
    _run_code(code, mod_globals, init_globals,
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "c:\Users\Damdae\OneDrive - SNU\Graduate\Project\Design\ITER\script\pyomo_test\ex2_simulator2.py", line 78, in <module>
    sim, tsim, profiles = simulate_model(model)
  File "c:\Users\Damdae\OneDrive - SNU\Graduate\Project\Design\ITER\script\pyomo_test\ex2_simulator2.py", line 45, in simulate_model
    tsim, profiles = sim.simulate(numpoints=100, integrator='vode')
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\pyomo\dae\simulator.py", line 891, in simulate
    tsim, profile = self._simulate_with_scipy(initcon, tsim, switchpts,
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\pyomo\dae\simulator.py", line 933, in _simulate_with_scipy
    profilestep = scipyint.integrate(tsim[i])
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\scipy\integrate\_ode.py", line 433, in integrate
    self._y, self.t = mth(self.f, self.jac or (lambda: None),
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\scipy\integrate\_ode.py", line 1009, in run
    y1, t, istate = self.runner(*args)
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\pyomo\dae\simulator.py", line 644, in _rhsfun
    residual.append(rhsdict[d]())
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\pyomo\core\base\param.py", line 829, in __call__
    return super(ScalarParam, self).__call__(exception=exception)
  File "C:\Users\Damdae\anaconda3\envs\arena\lib\site-packages\pyomo\core\base\param.py", line 168, in __call__
    raise ValueError(
ValueError: Error evaluating Param value ('eq1[{t}]'):
        The Param value is currently set to an invalid value.  This is
        typically from a scalar Param or mutable Indexed Param without
        an initial or default value.

Information on your system

Pyomo version: 6.23 Python version: 3.8.12 Operating system: Windows 10 How Pyomo was installed (PyPI, conda, source): pip Solver (if applicable): None (used Simulator)

Additional information

None

[blnicho]

@damdaepark thanks for reporting this. It's a known issue that I've been working on and I'm planning to have a fix for it in the next release.