Failure Of Basic Wall Constraint For Spring Mass Example

[SoundsSerious]

Description

I am evaluating dymos for evaluation of optimization of a dynamic hydraulic system operating in resonance.

I see the caveats regarding the system and variables must be continuous and differentiable so I though it would be a good idea to test dymos with something like an actuator limit for the spring mass example.

Adding a basic lower limit fails at x=-4. I am thinking whats going on is that the v paramater should be bounded as well but the solver cant easily assess that since only constant (not piecewise) boundaries are considered, there isn't an easy way to say, velocity should be zero at the wall.

I am highly interested in this library but It doesn't seem ready for generally bounded problems, which unfortunately constitute most dynamic systems.

Example

Code Example:

import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt

from dymos.examples.oscillator.oscillator_ode import OscillatorODE

# Instantiate an OpenMDAO Problem instance.
prob = om.Problem()

# We need an optimization driver.  To solve this simple problem ScipyOptimizerDriver will work.
prob.driver = om.ScipyOptimizeDriver()

# Instantiate a Dymos Trajectory and add it to the Problem model.
traj = dm.Trajectory()
prob.model.add_subsystem('traj', traj)

# Instantiate a Phase and add it to the Trajectory.
phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=10))
traj.add_phase('phase0', phase)

# Tell Dymos that the duration of the phase is bounded.
phase.set_time_options(fix_initial=True, fix_duration=True)

# Tell Dymos the states to be propagated using the given ODE.
phase.add_state('x', fix_initial=True, rate_source='v', targets=['x'], units='m')
phase.add_state('v', fix_initial=True, rate_source='v_dot', targets=['v'], units='m/s')

# The spring constant, damping coefficient, and mass are inputs to the system that
# are constant throughout the phase.
phase.add_parameter('k', units='N/m', targets=['k'])
phase.add_parameter('c', units='N*s/m', targets=['c'])
phase.add_parameter('m', units='kg', targets=['m'])

phase.add_path_constraint('x',lower=-4)

# Since we're using an optimization driver, an objective is required.  We'll minimize
# the final time in this case.
phase.add_objective('time', loc='final')

# Setup the OpenMDAO problem
prob.setup()

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

prob.set_val('traj.phase0.states:x', 10.0)
prob.set_val('traj.phase0.states:v', 0.0)

prob.set_val('traj.phase0.parameters:k', 1.0)
prob.set_val('traj.phase0.parameters:c', 0.5)
prob.set_val('traj.phase0.parameters:m', 1.0)

# Now we're using the optimization driver to iteratively run the model and vary the
# phase duration until the final y value is 0.
prob.run_driver()

# Perform an explicit simulation of our ODE from the initial conditions.
sim_out = traj.simulate(times_per_seg=50)

# Plot the state values obtained from the phase timeseries objects in the simulation output.
t_sol = prob.get_val('traj.phase0.timeseries.time')
t_sim = sim_out.get_val('traj.phase0.timeseries.time')

states = ['x', 'v']
fig, axes = plt.subplots(len(states), 1)
for i, state in enumerate(states):
    sol = axes[i].plot(t_sol, prob.get_val(f'traj.phase0.timeseries.{state}'), 'o')
    sim = axes[i].plot(t_sim, sim_out.get_val(f'traj.phase0.timeseries.{state}'), '-')
    axes[i].set_ylabel(state)
axes[-1].set_xlabel('time (s)')
fig.legend((sol[0], sim[0]), ('solution', 'simulation'), loc='lower right', ncol=2)
plt.tight_layout()
plt.show()

Dymos Version

1.10.0

Relevant environment information

No response

[robfalck]

Hello and thanks for checking out dymos.

The oscillator example posted is an example of solving an initial value problem in dymos. The variables that govern the trajectory are all fixed:

• The initial time and duration of the trajectory
• The initial values of the position x and velocity y
• The value of the parameters of the ODE: k, c, and m.

With all of these variables fixed, there is only a single physical path that the system can possibly take, so there are no degrees of freedom or "knobs" that can be turned to satisfy the wall path constraint at x=-4.

This example is a bit odd, because we're using an optimizer to satisfy the pseudospectral method defect constraints. To make the optimizer happy, we also have to give it some objective. This example used final time, but the optimizer can't do anything to actually change it since the phase initial time and duration are fixed. If you check out opt_report.html in the reports directory generated when you run this, you'll see that there are 60 equality constraints and 60 design variables. That's an NxN system and theres only a single solution.

Instead, we can modify the problem and ask dymos to solve "What is the largest initial displacement the block can have such that it doesn't violate the wall constraint?"

The changes needed in the code are as follows. First, we free up the initial value of the displacement, x.

phase.add_state('x', fix_initial=False, rate_source='v', targets=['x'], units='m')

Next, we need to change the objective to maximize the initial value of x. Otherwise, dymos is going to slam the duration of the phase to the minimum allowable value and not guarantee that the minimum displacement actually reaches x=-4.

phase.add_objective('x', loc='initial', scaler=-1.)

Also, I added the following line after the driver has been instantiated. This speeds up the calculation of the total derivatives needed by the optimizer:

prob.driver.declare_coloring()

Note that path constraints are only enforced at the discrete nodes. There may be some minor violation of the path constraint between nodes. You can get more nodes by using more segments or higher-order segments.

If you absolutely need to know the exact condition when the wall is reached, that would be a use case for a two phase problem and boundary constraints. The first phase would proceed for some minimum duration and have a final boundary constraint of v=0. With the appropriate problem formulation, you could also constrain that the final displacement x=-4.

If you're looking for a simulation tool that will "bounce" off of the walls, you could pose that in dymos with appropriate phases and boundary constraints but that's not really what it's designed to do. Phase-based specification of trajectories means you need some a-priori knowledge of the discontinuities in the states, such as the jump in velocity when the mass hits the wall.

import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt

from dymos.examples.oscillator.oscillator_ode import OscillatorODE

# Instantiate an OpenMDAO Problem instance.
prob = om.Problem()

# We need an optimization driver.  To solve this simple problem ScipyOptimizerDriver will work.
prob.driver = om.ScipyOptimizeDriver()
prob.driver.declare_coloring()

# Instantiate a Dymos Trajectory and add it to the Problem model.
traj = dm.Trajectory()
prob.model.add_subsystem('traj', traj)

# Instantiate a Phase and add it to the Trajectory.
phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=10))
traj.add_phase('phase0', phase)

# Tell Dymos that the duration of the phase is bounded.
phase.set_time_options(fix_initial=True, fix_duration=True)

# Tell Dymos the states to be propagated using the given ODE.
phase.add_state('x', fix_initial=False, rate_source='v', targets=['x'], units='m')
phase.add_state('v', fix_initial=True, rate_source='v_dot', targets=['v'], units='m/s')

# The spring constant, damping coefficient, and mass are inputs to the system that
# are constant throughout the phase.
phase.add_parameter('k', units='N/m', targets=['k'])
phase.add_parameter('c', units='N*s/m', targets=['c'])
phase.add_parameter('m', units='kg', targets=['m'])

phase.add_path_constraint('x',lower=-4)

# Since we're using an optimization driver, an objective is required.  We'll minimize
# the final time in this case.
phase.add_objective('x', loc='initial', scaler=-1.0)

# Setup the OpenMDAO problem
prob.setup()

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

prob.set_val('traj.phase0.states:x', 10.0)
prob.set_val('traj.phase0.states:v', 0.0)

prob.set_val('traj.phase0.parameters:k', 1.0)
prob.set_val('traj.phase0.parameters:c', 0.5)
prob.set_val('traj.phase0.parameters:m', 1.0)

# Now we're using the optimization driver to iteratively run the model and vary the
# phase duration until the final y value is 0.
prob.run_driver()

# Perform an explicit simulation of our ODE from the initial conditions.
sim_out = traj.simulate(times_per_seg=50)

# Plot the state values obtained from the phase timeseries objects in the simulation output.
t_sol = prob.get_val('traj.phase0.timeseries.time')
t_sim = sim_out.get_val('traj.phase0.timeseries.time')

states = ['x', 'v']
fig, axes = plt.subplots(len(states), 1)
for i, state in enumerate(states):
    sol = axes[i].plot(t_sol, prob.get_val(f'traj.phase0.timeseries.{state}'), 'o')
    sim = axes[i].plot(t_sim, sim_out.get_val(f'traj.phase0.timeseries.{state}'), '-')
    axes[i].set_ylabel(state)
axes[-1].set_xlabel('time (s)')
fig.legend((sol[0], sim[0]), ('solution', 'simulation'), loc='lower right', ncol=2)
plt.tight_layout()
plt.show()import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt

from dymos.examples.oscillator.oscillator_ode import OscillatorODE

# Instantiate an OpenMDAO Problem instance.
prob = om.Problem()

# We need an optimization driver.  To solve this simple problem ScipyOptimizerDriver will work.
prob.driver = om.ScipyOptimizeDriver()
prob.driver.declare_coloring()

# Instantiate a Dymos Trajectory and add it to the Problem model.
traj = dm.Trajectory()
prob.model.add_subsystem('traj', traj)

# Instantiate a Phase and add it to the Trajectory.
phase = dm.Phase(ode_class=OscillatorODE, transcription=dm.Radau(num_segments=10))
traj.add_phase('phase0', phase)

# Tell Dymos that the duration of the phase is bounded.
phase.set_time_options(fix_initial=True, fix_duration=True)

# Tell Dymos the states to be propagated using the given ODE.
phase.add_state('x', fix_initial=False, rate_source='v', targets=['x'], units='m')
phase.add_state('v', fix_initial=True, rate_source='v_dot', targets=['v'], units='m/s')

# The spring constant, damping coefficient, and mass are inputs to the system that
# are constant throughout the phase.
phase.add_parameter('k', units='N/m', targets=['k'])
phase.add_parameter('c', units='N*s/m', targets=['c'])
phase.add_parameter('m', units='kg', targets=['m'])

phase.add_path_constraint('x',lower=-4)

# Since we're using an optimization driver, an objective is required.  We'll minimize
# the final time in this case.
phase.add_objective('x', loc='initial', scaler=-1.0)

# Setup the OpenMDAO problem
prob.setup()

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

prob.set_val('traj.phase0.states:x', 10.0)
prob.set_val('traj.phase0.states:v', 0.0)

prob.set_val('traj.phase0.parameters:k', 1.0)
prob.set_val('traj.phase0.parameters:c', 0.5)
prob.set_val('traj.phase0.parameters:m', 1.0)

# Now we're using the optimization driver to iteratively run the model and vary the
# phase duration until the final y value is 0.
prob.run_driver()

# Perform an explicit simulation of our ODE from the initial conditions.
sim_out = traj.simulate(times_per_seg=50)

# Plot the state values obtained from the phase timeseries objects in the simulation output.
t_sol = prob.get_val('traj.phase0.timeseries.time')
t_sim = sim_out.get_val('traj.phase0.timeseries.time')

states = ['x', 'v']
fig, axes = plt.subplots(len(states), 1)
for i, state in enumerate(states):
    sol = axes[i].plot(t_sol, prob.get_val(f'traj.phase0.timeseries.{state}'), 'o')
    sim = axes[i].plot(t_sim, sim_out.get_val(f'traj.phase0.timeseries.{state}'), '-')
    axes[i].set_ylabel(state)
axes[-1].set_xlabel('time (s)')
fig.legend((sol[0], sim[0]), ('solution', 'simulation'), loc='lower right', ncol=2)
plt.tight_layout()
plt.show()

[SoundsSerious]

@robfalck Thanks for taking the time to work through this issue and explain some possible solutions to this problem to me.

I recognize dymos has a unique approach to dynamic problems, solving the trajectory as a state. That is very novel and I completely understand that this conventional dynamics problem isn't really a good fit for the library's spectral method. There will be plenty of scenarios where I expect the physics of my problem to rest on the actuator limit.

I understand there is an explicit shooting mode to dymos as well, I would be very interested in a way to use more complicated path boundaries, such as conditional wall boundaries / actuator limits in openmdao.

This would allow dymos to be used for general dynamics problems, in tandem with openmdao. Before I encountered dymos I was wondering how to combine openmdao with the time domain for optimization, and I recognize its a tricky intersection of solver / optimizer approaches, and one that I think dymos approaches intelligently given all the challenges of integrating the time domain into openmdao.

With that I'll leave it up to the maintainers of dymos to close this issue or leave it open in the case there might be a solution for more complex boundary conditions.

[robfalck]

There's a paper on the "sweeping gradient method" by Ben Margolis that accomplishes what you want, ODE integration with adjoint differentiation and discrete event detection. https://www.researchgate.net/publication/374475723_A_Sweeping_Gradient_Method_for_Ordinary_Differential_Equations_with_Events

I haven't implemented it in dymos yet. It's very much a "controls view" of solving the problem where things like path and boundary constraints are enforced by a controler modeled in the ODE,

[SoundsSerious]

Thats quite interesting! Thanks for sharing.

That would make dymos quite powerful with the ability to create an optimal path and controller! I think openmdao would probably be quite useful here as well with partial derivatives baked in.

I think thats probably quite a hard problem to generalize along with the existing methods of dymos.

I mentioned dymos in the python-controls library, as they had an open issue for a better system simulator / solver for integrated components so they might be interested in some kind of integration. Their library might assist with #1031