Use Dymos and OpenMDAO to simulate a pressurized bottle fluid expulsion

The primary issue that you have asked about, related to the No such ODE input error, is cased by the way you coded your ODE and more specifically the way you promoted variables and then added the ODE to the phase.

For example, you promoted your input P then set the state target to pressure_rate.P. This is incorrect. When you promoted P that effectively moved the name address up to the top level of the ODE, so the name target is just P now. You can read more about promotion vs connection in the docs. You have this problem in most of your script, where you are not accounting for promotion when you set targets.

Unfortunately, this is not the only issue in your script. There are several more, and enough that I am not able to get things fully working. Here are some other modest issues in rough order of importance:

1. The VolumeFlowRate component input and outputs are scalar, but seem to be intended to connect to the vector (of size num_nodes) variables of PressureRate. I suspect you meant to make them vector as well, but am not 100% sure
2. You have an execution order issue between PressureRate and VolumeRate. Pressure rate seems to need as an input Vl_dot, which comes from VolumeRate`, but you have added it first so it will run BEFORE the component providing its input value.
3. You had a typo in your set_val calls (explusion vs expulsion)
4. You did not have a deltap key in the parameter diction, but you did have a variable for it.

After fixing those, I could get the problem to start running but it did not converge or give an answer. You had solve_segments set to forward and had set 50 segments. Both of those seemed like bad settings to me, so I changed them to 3 segments, and removed the solve_segments option.

Then I was able to get the optimizer to take a few steps, but it errored with

            Current function value: [200.]
            Iterations: 6
            Function evaluations: 12
            Gradient evaluations: 2
Optimization FAILED.
Positive directional derivative for linesearch

Which indicated a problem with the derivatives. So I changed your setting for partial derivatives from fd to cs. That allowed it to iterate more, but still didn't converge. Without diving more into the physics of your problem I can't easily diagnose this further. I suspect you have some bad boundary conditions and probably bad initial guesses though.

Here is the modified script I came up with to at least get the optimizer iterating.

import numpy as np
import openmdao.api as om 
import dymos as dm

class PressureRate(om.ExplicitComponent):

    def initialize(self):
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']

        # Inputs
        self.add_input('p', shape=(nn,), desc='Pressure inside the nox bottle', units='Pa')
        self.add_input('Vb', shape=(nn,), desc='Bottle volume', units='m**3')
        self.add_input('Vl', shape=(nn,), desc='Liquid volume', units='m**3')
        self.add_input('Vl_dot', shape=(nn,), desc='Liquid volume flow rate', units='m**3/s')
        self.add_input('gamma', shape=(nn,), desc='Heat capacity ratio')

        # Outputs
        self.add_output('p_dot', val=np.ones(nn), desc='Pressure change rate', units='Pa/s')
        self.declare_partials(of='*', wrt='*', method='cs')

    def compute(self, inputs, outputs):
        p = inputs['p']
        Vb = inputs['Vb']
        Vl = inputs['Vl']
        Vl_dot = inputs['Vl_dot']
        gamma = inputs['gamma']

        outputs['p_dot'] = gamma * p/(Vb - Vl) * Vl_dot

class VolumeFlowRate(om.ExplicitComponent):
    """
    A Dymos ODE for a damped harmonic oscillator.
    """

    def initialize(self):
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']

        # Inputs
        self.add_input('p', shape=(nn,), desc='Pressure inside the nox_bottle', units='Pa')
        self.add_input('pout', shape=(nn,), desc='Pressure outside the nox_bottle', units='Pa')
        self.add_input('deltap', shape=(nn,), desc='Nox bottle pressure losses', units='Pa')
        self.add_input('rhol', shape=(nn,), desc='Liquid density', units='kg/m**3')
        self.add_input('Aout', shape=(nn,), desc='Output nox_bottle area', units='m**2')

        # Outputs
        self.add_output('Vl_dot', shape=(nn,), desc='Volume flow rate', units='m**3/s')

        self.declare_partials(of='*', wrt='*', method='cs')


    def compute(self, inputs, outputs):
        p = inputs['p']
        pout = inputs['pout']
        deltap = inputs['deltap']
        rhol = inputs['rhol']
        Aout = inputs['Aout']

        outputs['Vl_dot'] = Aout*np.sqrt(2/rhol*(p - pout - deltap))

class BottleModelODE(om.Group):

    def initialize(self):
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']

        self.add_subsystem('volume_flow_rate', subsys=VolumeFlowRate(num_nodes=nn),
                           promotes_inputs=['p', "pout", 'deltap', 'rhol', "Aout"], promotes_outputs=['Vl_dot'])
        self.add_subsystem('pressure_rate', subsys=PressureRate(num_nodes=nn),
                           promotes_inputs=['p', "Vb", "Vl", "Vl_dot", "gamma"], promotes_outputs=['p_dot'])
        



def expulsion_phase_fn(transcription: dm.transcriptions.pseudospectral.radau_pseudospectral.Radau, pamb: float):

    phase = dm.Phase(ode_class=BottleModelODE, transcription=transcription)

    phase.set_time_options(fix_initial=True, fix_duration=True)

    phase.add_state('p', units='bar', rate_source='p_dot',
                    targets=['p'], fix_initial=True, fix_final=False, lower=pamb)
    phase.add_state('Vl', units='m**3', rate_source='Vl_dot', targets=['Vl'],
                    fix_initial=True, fix_final=False, lower=0)

    phase.add_parameter('Vb', targets=['Vb'], units='m**3')
    phase.add_parameter('gamma', targets=['gamma'])
    phase.add_parameter('rhol', targets=['rhol'], units='kg/m**3')
    phase.add_parameter('Aout', targets=['Aout'], units='m**2')
    phase.add_parameter('pout', targets=['pout'], units="Pa")
    phase.add_parameter('deltap', targets=['deltap'], units="Pa")

    return phase

def trajectory(pamb: float):

    # transcript = dm.Radau(num_segments=50, solve_segments='forward')
    transcript = dm.Radau(num_segments=3)
    traj = dm.Trajectory()

    # Add phases to trajectory
    expulsion_phase = traj.add_phase('expulsion', expulsion_phase_fn(transcription=transcript, pamb=pamb))

    return traj, expulsion_phase



if __name__ == "__main__": 

    # Set ambiant conditions
    Tamb = 20 + 273.15
    pamb = 100*10**3
    deltap = 0
    Vb = 5*10**-3
    Aout = 10*10**-4

    # Set NOX bottle properties up
    bottle_params = {"Vb": 5*10**-3, "gamma": 1.4, "Aout": 3*10**-2, "rhol": 1000, "pout": 100*10**3, "pinit": 300*10**3, "Vl": 1*10**-3}

    # Instantiate an OpenMDAO Problem instance
    prob = om.Problem(model=om.Group())
    prob.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')

    # Instantiate a Dymos trjectory and add it to the Problem model
    traj, phase = trajectory(pamb=100*10*3)
    phase.add_objective("time", loc="final")

    # Setup the OpenMDAO problem
    prob.model.add_subsystem("traj", traj)
    prob.setup()

    # Assign values to the times and states
    prob.set_val('traj.expulsion.t_initial', 0.0)
    prob.set_val('traj.expulsion.t_duration', 200.0)

    prob.set_val('traj.expulsion.states:p', bottle_params["pinit"])
    prob.set_val('traj.expulsion.states:Vl', bottle_params["Vl"])

    prob.set_val('traj.expulsion.parameters:Vb', bottle_params["Vb"])
    prob.set_val('traj.expulsion.parameters:gamma', bottle_params["gamma"])
    prob.set_val('traj.expulsion.parameters:rhol', bottle_params["rhol"])
    prob.set_val('traj.expulsion.parameters:Aout', bottle_params["Aout"])
    prob.set_val('traj.expulsion.parameters:pout', bottle_params["pout"])
    prob.set_val('traj.expulsion.parameters:deltap', deltap)

    prob.run_driver()