使用Dymos和OpenMDAO模拟加压瓶流体排出

Question

我正在打开这个新线程，因为我正在寻找一些使用Dymos，以模拟一个动态系统。

事实上，我正在尝试模拟一个由一个加压的瓶子和里面的液体组成的系统。当t=0时，压力推动流体通过瓶出口，因此瓶内压力减小。我的目标是模拟瓶内压力的行为以及从瓶中逸出的流体体积流。我发现了一个Dymos的例子，它与我正在尝试做的非常相似，但更简单。https://openmdao.github.io/dymos/examples/water_rocket/water_rocket.html

为了模拟我的系统，我使用了两个明确的组件:PressureRate, VolumeFLowRate。然后我定义组组件PressureModelODE来连接这最后两个组件和它们的变量。

以下是这些组件:

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='fd')
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']
def setup(self):
# Inputs
self.add_input('p', desc='Pressure inside the nox_bottle', units='Pa')
self.add_input('pout', desc='Pressure outside the nox_bottle', units='Pa')
self.add_input('deltap', desc='Nox bottle pressure losses', units='Pa')
self.add_input('rhol', desc='Liquid density', units='kg/m**3')
self.add_input('Aout', desc='Output nox_bottle area', units='m**2')
# Outputs
self.add_output('Vl_dot', desc='Volume flow rate', units='m**3/s')
self.declare_partials(of='*', wrt='*', method='fd')
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('pressure_rate', subsys=PressureRate(num_nodes=nn),
promotes_inputs=['p', "Vb", "Vl", "Vl_dot", "gamma"], promotes_outputs=['p_dot'])
self.add_subsystem('volume_flow_rate', subsys=VolumeFlowRate(num_nodes=nn),
promotes_inputs=['p', "pout", 'deltap', 'rhol', "Aout"], promotes_outputs=['Vl_dot'])
self.connect('pressure_rate.p', 'volume_flow_rate.p')
self.connect('pressure_rate.Vl_dot', 'volume_flow_rate.Vl_dot')

然后，为了求解这些方程并模拟我的模型，我基于振荡器的例子构建了一个程序:https://openmdao.github.io/dymos/getting_started/intro_to_dymos/intro_segments.html

我正在定义一个叫做"爆炸"的阶段。通过使用以下函数:

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='pressure_rate.p_dot',
targets=['pressure_rate.p', "volume_flow_rate.p"], fix_initial=True, fix_final=False, lower=pamb)
phase.add_state('Vl', units='m**3', rate_source='volume_flow_rate.Vl_dot', targets=['pressure_rate.Vl'],
fix_initial=True, fix_final=False, lower=0)
phase.add_parameter('Vb', targets=['pressure_rate.Vb'], units='m**3')
phase.add_parameter('gamma', targets=['pressure_rate.gamma'])
phase.add_parameter('rhol', targets=['volume_flow_rate.rhol'], units='kg/m**3')
phase.add_parameter('Aout', targets=['volume_flow_rate.Aout'], units='m**2')
phase.add_parameter('pout', targets=['volume_flow_rate.pout'], units="Pa")
phase.add_parameter('deltap', targets=['volume_flow_rate.deltap'], units="Pa")
return phase

然后，我用这个函数定义一个轨迹:

def trajectory(pamb: float):
transcript = dm.Radau(num_segments=50, solve_segments='forward')
traj = dm.Trajectory()
# Add phases to trajectory
expulsion_phase = traj.add_phase('expulsion', 
expulsion_phase_fn(transcription=transcript, pamb=pamb))
return traj, expulsion_phase

最后，我正在设置OpenMDAO问题，提供初始值，…通过执行以下基于振荡器示例的行:

def launch_compt():
# 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.explusion.t_initial', 0.0)
prob.set_val('traj.explusion.t_duration', 200.0)
prob.set_val('traj.explusion.states:p', bottle_params["pinit"])
prob.set_val('traj.explusion.states:Vl', bottle_params["Vl"])
prob.set_val('traj.explusion.parameters:Vb', bottle_params["Vb"])
prob.set_val('traj.explusion.parameters:gamma', bottle_params["gamma"])
prob.set_val('traj.explusion.parameters:rhol', bottle_params["rhol"])
prob.set_val('traj.explusion.parameters:Aout', bottle_params["Aout"])
prob.set_val('traj.explusion.parameters:pout', bottle_params["pout"])
prob.set_val('traj.explusion.parameters:deltap', bottle_params["deltap"])
prob.run_driver()

不幸的是，这不起作用，我不明白为什么。它返回给我，参数Vb(瓶总量)没有提供，但我不明白为什么:它是在我添加参数到问题时提供的，就像在振荡器示例中一样。

在这方面，我正在联系，希望能找到一些帮助。谢谢你的回答。

PS:这是我得到的错误信息，当我试图运行程序:

raise ValueError(f'Invalid parameter in phase `{self.pathname}`.n{str(e)}') from e
ValueError: Invalid parameter in phase `traj.phases.expulsion`.
Parameter `Vb` has invalid target(s).
No such ODE input: 'pressure_rate.Vb'.

Answer 1

您所询问的与No such ODE input错误相关的主要问题是您编码ODE的方式，更具体地说，是您提升变量并将ODE添加到阶段的方式。

例如，您将输入提升为P，然后将状态目标设置为pressure_rate.P。这是不正确的。当您提升P时，它有效地将名称地址移动到ODE的顶层，因此名称目标现在只是P。你可以在文档中阅读更多关于推广与联系的内容。你的大部分脚本都有这个问题，当你设定目标时，你没有考虑到促销。

不幸的是，这不是脚本中唯一的问题。还有几个，而且足够多，我不能让事情完全工作。以下是其他一些适度的问题，按其重要性大致排序:

1. VolumeFlowRate组件的输入和输出是标量，但似乎打算连接到PressureRate的向量(大小为num_nodes)变量。我怀疑你的意思是让他们矢量，但我不是100%确定
2. PressureRate和VolumeRate之间有执行顺序问题。Pressure rate seems to need as an inputVl_dot, which comes fromVolumeRate '，但是你已经添加了它，所以它将在组件提供其输入值之前运行。
3. 你在set_val调用中有一个错字(explusionvsexpulsion)
4. 您的参数中没有deltap键，但是您有一个变量。

修复这些后，我可以让问题开始运行，但它没有收敛或给出答案。您将solve_segments设置为转发，并设置了50个段。这两种设置对我来说都很糟糕，所以我把它们改成了3段，并删除了solve_segments选项。

然后我可以让优化器执行几个步骤，但是

出错了。

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

这表明衍生品有问题。所以我把偏导数的设置从fd改成了cs。这使得它可以迭代更多，但仍然没有收敛。如果不深入研究你的问题的物理原理，我就无法轻易地进一步诊断。我怀疑你有一些不好的边界条件，可能是错误的初步猜测。

下面是我修改后的脚本，至少可以让优化器迭代。

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()