Code Monkey home page Code Monkey logo

Comments (4)

SoundsSerious avatar SoundsSerious commented on August 26, 2024 1

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

from dymos.

robfalck avatar robfalck commented on August 26, 2024

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

from dymos.

SoundsSerious avatar SoundsSerious commented on August 26, 2024

@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.

from dymos.

robfalck avatar robfalck commented on August 26, 2024

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,

from dymos.

Related Issues (20)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.