Water rocket example

[bbahiam]

Summary of Issue

Hi everyone,

I am trying to create a water rocket example for dymos by modifying the cannonball example but I am running into some problems (mostly due to my inexperience with dymos). I would really appreciate a hand!

Issue Type

• [ ] Bug
• [x] Enhancement

Description

My idea to create the water rocket example is to add an extra propelled phase before the cannonball ascent phase. This phase has two extra states compared to cannonball: the pressure and the water volume inside the rocket. It runs until the water inside the rocket is depleted.

So far I have created a group for the "water engine" and added it to a modified version of the CannonballODE. Here is a N2 diagram of it: n2.zip

Now I am trying to link it with the other phases, but I am running into a error, as explained below.

Example

The modified example is in this file: https://github.com/bbahiam/dymos/blob/water_rocket_example/dymos/examples/water_rocket/doc/test_doc_water_rocket.py

I am running into the following error

Traceback (most recent call last):
  File "doc/test_doc_water_rocket.py", line 141, in test_water_rocket_for_docs
    p.setup()
  File "OpenMDAO/openmdao/core/problem.py", line 758, in setup
    derivatives, self.options)
  File "OpenMDAO/openmdao/core/system.py", line 790, in _setup
    self._setup_procs(self.pathname, comm, mode, self._problem_options)
  File "OpenMDAO/openmdao/core/group.py", line 420, in _setup_procs
    subsys._setup_procs(subsys.pathname, sub_comm, mode, prob_options)
  File "OpenMDAO/openmdao/core/group.py", line 361, in _setup_procs
    self.setup()
  File "dymos/dymos/trajectory/trajectory.py", line 417, in setup
    self._setup_design_parameters()
  File dymos/dymos/trajectory/trajectory.py", line 311, in _setup_design_parameters
    if tgts[phase_name] is None:
KeyError: 'propelled_ascent'

Can you give me any pointers on what I am missing here?

Thanks!

Environment

Operating System: Ubuntu 18.04.1 Python environment: Python 3.6.9 Packages:

-e git+git@github.com:bbahiam/dymos.git@f6a86c105b7b14c038dff1c36a7bbb74ac7ddfcb#egg=dymos
-e git+git@github.com:OpenMDAO/OpenMDAO.git@073f1a86171769eff7e9a7bc380c88fa11c15611#egg=openmdao
-e git+git@github.com:mdolab/pyoptsparse.git@89496fcc78476f497332184fd12d4d001f6c6b9c#egg=pyoptsparse

[JustinSGray]

could you post a gist of your run script, or a link to the specific commit of a repo?

[bbahiam]

I am working on a fork of dymos. This is the current commit https://github.com/bbahiam/dymos/commit/f6a86c105b7b14c038dff1c36a7bbb74ac7ddfcb

Everything I did was in dymos/examples/water_rocket and the run script is in dymos/examples/water_rocket/doc/test_doc_water_rocket.py.

[robfalck]

So this is exposing a few issues in the parameter handling for trajectories in Dymos. In the meantime, I suggest the following:

1. Provide targets for all trajectory input and design parameters

        # Add internally-managed design parameters to the trajectory.
        traj.add_design_parameter('CD',
                                  targets={'propelled_ascent': ['aero.CD'],
                                           'ballistic_ascent': ['aero.CD'],
                                           'descent': ['aero.CD'],},
                                  val=0.5, units=None, opt=False)
        traj.add_design_parameter('CL',
                                  targets={'propelled_ascent': ['aero.CL'],
                                           'ballistic_ascent': ['aero.CL'],
                                           'descent': ['aero.CL'],},
                                  val=0.0, units=None, opt=False)
        traj.add_design_parameter('alpha',
                                  targets={'propelled_ascent': ['eom.alpha'],
                                           'ballistic_ascent': ['eom.alpha'],
                                           'descent': ['eom.alpha'],},
                                  val=0.0, units='deg', opt=False)

and

        traj.add_input_parameter('m', units='kg', val=1.0,
                                 targets={'propelled_ascent': None,
                                          'ballistic_ascent': 'mass',
                                          'descent': 'mass'})

We'll fix dymos to be a bit better about handling missing phases from the targets list and some better debugging messages.

Having fixed those, there are next a few issues in the ODE.

2. In water_engine_comp.py there you promote the names Vdot but then issue the explicit connection

self.connect('water_flow_rate.Vdot', 'pressure_rate.Vdot')

Since these variables have been promoted, they will implicitly be connected and this connect statement shouldn't be issued.

Theres also an issue with a connection inside the cannonball ode (kinetic_energy.v). The cannonball example is working on dymos/master but not on your fork, so I'm not sure what's changed here.

[robfalck]

Connection issues pointed out here should be addressed by #325

[bbahiam]

Thanks for the pointers Rob. I am working on it now and will let you guys know what I can or can't figure out

[bbahiam]

Hi, I managed to evolve a bit. I merged Rob's pull and I think I figured out all connection issues. Now the optimization is running but it is diverging. I will try to isolate the propelled phase to see what is going on.

The current commit is https://github.com/bbahiam/dymos/commit/ad56d5e3e020d3f11aa739a64f17340a78636b14

[robfalck]

Thanks, I'm going to take a look at this when I get a chance.

[bbahiam]

No rush, I am making good progress on my own here. I isolated the phase and the optimization converges now. The next step is to plot the results and see how they look. I will let you know if I need more help.

Thanks

[bbahiam]

Is there any way to simulate the ODE without optimizing it? I have tried running traj.simulate() before solving the problem, but I get an error:

RuntimeError: Trajectory (traj): Unable to list outputs on a Group until model has been run.

I also could not find the documentation for this function:

(Pdb) help(traj.simulate)
*** No help for '(traj.simulate)`

EDIT: I found the documentation at https://openmdao.github.io/dymos/feature_reference/phases/simulate/simulate.html

[robfalck]

Call run_model() (rather than run_driver()) and it will just propagate the calculations on the system, and then simulate will be able to run.

[bbahiam]

Thanks, it works.

[robfalck]

Any plans to share this model?

[bbahiam]

I was planning to push it to the main dymos repository once it is working and documented as well as the other examples, if you are happy with that.

On Fri, May 1, 2020, 15:19 Rob Falck notifications@github.com wrote:

Any plans to share this model?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/OpenMDAO/dymos/issues/322#issuecomment-622525498, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANOXUPOZVBXYZQN2XNTX5ELRPMOFLANCNFSM4MTKGG3Q .

[robfalck]

Thanks, would be great. It also makes for a pretty accessible example of optimal control for STEM.

[bbahiam]

Yeah, actually the idea for this came about because I remembered that in high school I did a project on optimizing the amount of water in a water rocket, so I thought why not do it again but computationally.

The final goal is to optimize the water volume in the bottle and the empty weight.

On Fri, May 1, 2020, 16:09 Rob Falck notifications@github.com wrote:

Thanks, would be great. It also makes for a pretty accessible example of optimal control for STEM.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/OpenMDAO/dymos/issues/322#issuecomment-622544990, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANOXUPIP2JBTFLQHUBNIPSLRPMT5RANCNFSM4MTKGG3Q .

[bbahiam]

Happy to say that both the isolated propelled ascent phase and the full ascent optimizing for height are working. This is in commit e3f6622 "Height optimization working, but ill-conditioned". Here are some plots

Internal pressure:

Water volume

Velocity (propelled phase)

Height (propelled phase)

The problem that I am hitting now is that the 2D EOM are ill-conditioned at the final point of the optimum solution for the height problem, since the velocity is zero hence the flight path angle is not defined. Since I only need 1D dynamics for this problem, I am thinking about removing some states (gam and r). Is that possible? Can I just set them to a constant value?

[robfalck]

For the sake of this problem, you could just start the rocket off with some slightly off-vertical flight path angle. I believe theres a singularity in the range rate derivatives at 90 deg flight path angle. If started off-vertical the velocity should never completely reach zero.

Better yet is to use cartesian coordinate equations of motion, which don't suffer from the singularity at v=0 or gam=90. The single stage to orbit example already uses these: https://github.com/OpenMDAO/dymos/blob/master/dymos/examples/ssto/launch_vehicle_2d_eom_comp.py

[bbahiam]

Setting the Flight path angle close to vertical is working well enough for now. I am getting nice optimization results using range and height as objective functions:

	range	height
maximum range	97.0779 m	0.6767 m
maximum height	26.0491 m	60.9630 m
maximum velocity	43.4167 m/s	47.5811 m/s
launch angle	36.8861 deg	89.5302 deg
empty mass	0.1681 kg	0.1262 kg
water volume	1.0149 L	1.0149 L

data obtained from https://github.com/bbahiam/dymos/commit/a56e4123806fd8c4a33fbbd26582ba3461436fbd

Now, is there any way to do a DOE on Dymos? I want to vary the empty mass and initial water volume while satisfying the ODE. I was thinking about using OpenMDAO's DOEDriver, but I think that will not enforce the ODE collocation constraints. Is that right? Any suggestions? The end goal is a contour plot of the maximum height with empty mass and initial water volume as the axis.

[robfalck]

Nice work.

Yeah in thinking about it more, flight path angle causes a singularity in azimuth which isn't included in these simple 2D EOM.

I personally haven't used the DOEDriver. For a basic approach I'd wrap a loop around the driver execution that set_val on the appropriate variables before each solve.

[bbahiam]

Yeah, I though the singularity was the case when the optimizer didn't want to go more vertical than 89.5deg.

Cool, I was thinking about writing that loop, but just wanted to make sure I was not missing a cleaner way

Cheers

On Wed, May 6, 2020, 11:12 Rob Falck notifications@github.com wrote:

Nice work.

Yeah in thinking about it more, flight path angle causes a singularity in azimuth which isn't included in these simple 2D EOM.

I personally haven't used the DOEDriver. For a basic approach I'd wrap a loop around the driver execution that set_val on the appropriate variables before each solve.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/OpenMDAO/dymos/issues/322#issuecomment-624708774, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANOXUPNGCIZRGG5BPSVJUW3RQF455ANCNFSM4MTKGG3Q .

[robfalck]

The fact that you only experience problems as velocity approached zero also indicates that the 90 deg flight path angle is not itself the cause of the singularity, but zero velocity is. In this EOM theres no such thing as negative velocity so it's discontinuous in first derivative at v=0.

[bbahiam]

Can I by-pass the optimizer somehow and still satisfy the constraints generated by the ODE? I tried these changes

-        p.model.linear_solver = om.DirectSolver()
+        p.model.non_linear_solver = om.NewtonSolver()
...
-            p.run_problem()
+            p.run_model()

But it did not work, i.e. the ODE was not satisfied. Basically I want to have a newton solver responsible for the collocation constraints so that I can run a MDA. Of course this is the same as using the optimizer with enough constraints so that the feasible region is a single point, but I would like to know if the dymos implementation supports doing that in a more explicit way.

In the code, how are the constraints passed to the optimizer?

[robfalck]

So Dymos follows a lot of other pseudospectral software in that, by default, the constraints that enforce physical correctness (the defects) are imposed by constraints for the optimizer, along with the design variables (the values of the states and controls throughout the trajectory, and the static design parameters). This is why you need to use run_driver, since a single execution of run_model will just compute the defects but not drive them to zero.

BUT WAIT...

OpenMDAO makes it somewhat simple for us to treat the state history values as implicit variables instead of design variables. The corresponding residuals are just the defects. If you run Dymos in this mode, a simple run model will basically give you the time history that results from your assumed control profile, design parameters, and input parameters. To enable this, you can use the keyword argument solve_segments=True when you provide add_state or set_state_options. Or to do it for all states in the phase, provide solve_segments=True as an argument to your transcription.

This feature hasn't been documented because it's been experimental, but I'm confident enough that it will work on your problem. You can see how we setup the brachistochrone to optionally use solve_segments here: https://github.com/OpenMDAO/dymos/blob/master/dymos/examples/brachistochrone/test/ex_brachistochrone.py

[bbahiam]

But this still leaves the linking between the phases to the optimizer, right? i.e., it solves the collocation defects in each phase but not the linking defects. It works for the brachistochrone problem because it only has one phase... Or have I gotten it wrong?

In any case, this is not a big issue because I can work using the optimizer to solve the defects, I just set opt=False in all my design variables.

[robfalck]

You can inform the link_phases method that phases should be linked via connection, rather than by constraint: https://github.com/OpenMDAO/dymos/blob/master/dymos/examples/battery_multibranch/test/test_multibranch_trajectory.py#L166

[robfalck]

@bbahiam How is this going. Do you consider this issue to be closed for now? Will you be submitting an associated pull request for an example case?

[bbahiam]

Hey Rob, this is almost done, I am just working on the documentation now. I have refactored the code and barring some minor adjustments that part is done. This is my first time using sphinx too. I am trying to follow the code pattern of creating a .rst file under docs/examples and linking it to the code. You can close this for now and if I need more assistent I can get in touch again. I will submit a pull request, I have just been a bit short on time these last few weeks.

Cheers, Benardo

On Fri, Jun 26, 2020 at 10:50 AM Rob Falck notifications@github.com wrote:

@bbahiam https://github.com/bbahiam How is this going. Do you consider this issue to be closed for now? Will you be submitting an associated pull request for an example case?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/OpenMDAO/dymos/issues/322#issuecomment-650220492, or unsubscribe https://github.com/notifications/unsubscribe-auth/ANOXUPKILBSEEMRHP6NMYSTRYSYSRANCNFSM4MTKGG3Q .

[robfalck]

Sounds good Bernardo. Thank you for working through this.