array parameter not passing correctly in simulate()

akrdunn commented 3 years ago

Summary of Issue

A 1D array declared as a phase parameter with dynamic=False does not get passed correctly when calling traj.simulate()

Issue Type

[x] Bug
[ ] Enhancement
[ ] Docs
[ ] Miscellaneous

Description

I have an array that I need to pass into my ODE. I'm declaring the array with phase.add_parameter and dynamic=False. When calling dm.run_problem, I can confirm that the connected data gets passed in correctly and everything runs. When calling traj.simulate() however, it appears that the array gets connected improperly and the first value in the array is passed to every element in the ODE.

Example

I modified the brachistochrone example slightly to illustrate the issue: https://github.com/akrdunn/dymos_bug The print statements show the passed value of "array" at each function call. It's assigned correctly as [ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.] initially but changes to all ones in simulate.

Environment

Operating System: Windows 10 Python environment: Anaconda Python 3.8.3 Packages: Dymos 0.16.1, OpenMDAO 3.4.0

akrdunn commented 3 years ago

This issue is happening in phase.initialize_values_from_phase. Because I'm using openmdao 3.4.0 it's using val = phs.get_val(f'parameters:{name}')[0, ...] instead of just val = phs.get_val(f'parameters:{name}'). Removing the if statement and using the latter of those two fixes the problem.

robfalck commented 3 years ago

I just got around to testing this issue, and it appears that it already works correctly on the current master branch. (0.17.1-dev). The following test passes:

import unittest
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal
from openmdao.utils.testing_utils import use_tempdirs
import dymos as dm
from dymos.examples.brachistochrone import BrachistochroneODE
import numpy as np

@use_tempdirs
class TestSimulateArrayParam(unittest.TestCase):

    def test_simulate_array_param(self):
        #
        # Initialize the Problem and the optimization driver
        #
        p = om.Problem(model=om.Group())
        p.driver = om.ScipyOptimizeDriver()
        p.driver.declare_coloring()

        #
        # Create a trajectory and add a phase to it
        #
        traj = p.model.add_subsystem('traj', dm.Trajectory())

        phase = traj.add_phase('phase0',
                               dm.Phase(ode_class=BrachistochroneODE,
                                        transcription=dm.GaussLobatto(num_segments=10)))

        #
        # Set the variables
        #
        phase.set_time_options(initial_bounds=(0, 0), duration_bounds=(.5, 10))

        phase.add_state('x', fix_initial=True, fix_final=True, rate_source='xdot')

        phase.add_state('y', fix_initial=True, fix_final=True, rate_source='ydot')

        phase.add_state('v', fix_initial=True, fix_final=False, rate_source='vdot')

        phase.add_control('theta', continuity=True, rate_continuity=True,
                          units='deg', lower=0.01, upper=179.9)

        phase.add_parameter('g', units='m/s**2', val=9.80665)
        phase.add_parameter('array', units=None, shape=(10,), dynamic=False)

        # dummy array of data
        indeps = p.model.add_subsystem('indeps', om.IndepVarComp(), promotes=['*'])
        indeps.add_output('array', np.linspace(1, 10, 10), units=None)
        # add dummy array as a parameter and connect it
        p.model.connect('array', 'traj.phase0.parameters:array')

        #
        # Minimize time at the end of the phase
        #
        phase.add_objective('time', loc='final', scaler=10)

        p.model.linear_solver = om.DirectSolver()

        #
        # Setup the Problem
        #
        p.setup()

        #
        # Set the initial values
        #
        p['traj.phase0.t_initial'] = 0.0
        p['traj.phase0.t_duration'] = 2.0

        p['traj.phase0.states:x'] = phase.interpolate(ys=[0, 10], nodes='state_input')
        p['traj.phase0.states:y'] = phase.interpolate(ys=[10, 5], nodes='state_input')
        p['traj.phase0.states:v'] = phase.interpolate(ys=[0, 9.9], nodes='state_input')
        p['traj.phase0.controls:theta'] = phase.interpolate(ys=[5, 100.5], nodes='control_input')

        #
        # Solve for the optimal trajectory
        #
        dm.run_problem(p, simulate=True)

        # Test the results
        sol_results = om.CaseReader('dymos_solution.db').get_case('final')
        sim_results = om.CaseReader('dymos_solution.db').get_case('final')

        sol = sol_results.get_val('traj.phase0.timeseries.parameters:array')
        sim = sim_results.get_val('traj.phase0.timeseries.parameters:array')

        assert_near_equal(sol-sim, np.zeros_like(sol))

Note: a pending dymos update will allow parameters to be accessed from the save recorded files as traj.phase0.parameters:array instead of needing to pull them from the timeseries.

robfalck commented 3 years ago

@akrdunn if you could please test this with the latest master version (0.17.1-dev) I think it's already been resolved.

OpenMDAO / dymos