Comments (4)
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.
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 velocityy
- The value of the parameters of the ODE:
k
,c
, andm
.
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.
@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.
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)
- There is a bug in PETSc 3.21.1 that is affecting the Latest Test workflow
- Drop support for Python 3.8
- The matplotlib.cm.get_cmap() function has been removed as of mpl 3.9.0
- Unit ambiguity issue in Birkhoff transcription HOT 1
- Issue on page /examples/min_time_climb/min_time_climb.html HOT 1
- Potential Bug in Supersonic Interceptor MTTC Example HOT 1
- Dymos introspection fails if ODEs use OpenMDAO's runtime shape detection.
- input_initial breaks approximated total jacobeans
- The `in1d` function has been deprecated for NumPy 2.0. Replace with `isin`.
- Computational performance is significantly worse when using Birkhoff HOT 3
- Issue with ".set_time_val/.set_state_val/.set_control_val" in the Brachistochrone problem HOT 1
- `phase.load_case` should use the new `phase.set_xxx_vals` API
- Update tests on reports to match OpenMDAO API changes
- Change Birkhoff transcription interface to be more consistent with others.
- Unexpected behavior when OPENMDAO_REPORTS=0 HOT 1
- TimeseriesOutputComp allocates memory unnecessarily
- Bokeh timeseries plot hover shoudl show point value, not cursor location.
- Avoid use of `copy_build_artifacts.py` to copy over all outputs during doc build process.
- Remove use of `newshape` argument in numpy reshape for numpy 2.1 compatibility.
- Undefined function reference
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
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.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from dymos.