Search code examples
optimizationcontrolsopenmdaotranscriptionoptimal

How to specify integration of a quantity as an objective in Dymos


I am implementing Bryson-Denham problem. The objective is:

$J=\frac{1}{2}\int_{0}^{1}u^2\left(t\right)dt$

and in the doc of Dymos, all explanation and examples state objective value as a scalar at loc=initial or loc=final. I could not find any example that use integral of some quantity over time as an objective function. Is this possible? How can I implement this?

FYI, Bryson-Denham problem is well-explained in this page: https://www.gpops2.com/Examples/Bryson-Denham.html


Solution

  • Dymos will integrate any state you give it. In this case, you need to add a state for J and then also compute a state rate for it --- J_dot.

    import openmdao.api as om 
    import dymos as dm
    
    class BrysonDedhamODE(om.ExplicitComponent):
    
    
        def initialize(self):
            self.options.declare('num_nodes', types=int)
    
        def setup(self):
            nn = self.options['num_nodes']
    
            # static parameters
            self.add_input('x', shape=nn)
            self.add_input('v', shape=nn)
            self.add_input('u', shape=nn)
            self.add_input('J', shape=nn)
    
    
    
            # state rates
            self.add_output('x_dot', shape=nn, tags=['dymos.state_rate_source:x'])
            self.add_output('v_dot', shape=nn, tags=['dymos.state_rate_source:v'])
            self.add_output('J_dot', shape=nn, tags=['dymos.state_rate_source:J'])
    
            # Ask OpenMDAO to compute the partial derivatives using complex-step
            # with a partial coloring algorithm for improved performance, and use
            # a graph coloring algorithm to automatically detect the sparsity pattern.
            self.declare_coloring(wrt='*', method='cs')
    
    
    
        def compute(self, inputs, outputs):
    
            v, u, j = inputs["v"], inputs["u"], inputs["J"]
    
            outputs['x_dot'] = v        
            outputs['v_dot'] = u        
            outputs['J_dot'] = 0.5*u**2
    
    
    p = om.Problem()
    
    p.driver = om.pyOptSparseDriver()
    p.driver.options['optimizer'] = 'SLSQP'
    p.driver.declare_coloring()
    
    traj = p.model.add_subsystem('traj', dm.Trajectory())
    
    transcription = dm.Radau(num_segments=10, order=3, compressed=True)
    phase0 = dm.Phase(ode_class=BrysonDedhamODE, transcription=transcription)
    
    traj.add_phase('phase0', phase0)
    
    
    phase0.set_time_options(fix_initial=True, fix_duration=True)
    
    phase0.set_state_options("x", fix_initial=True, fix_final=True, lower=0, upper=2)
    phase0.set_state_options("v", fix_initial=True, fix_final=True, lower=-2, upper=2)
    phase0.set_state_options("J", fix_initial=False, fix_final=False,lower=-10, upper=10)
    
    phase0.add_control('u', lower=-10, upper=10)
    phase0.add_path_constraint('x', upper=1/9.)
    
    phase0.add_objective('J', loc="final")
    
    p.setup()
    
    #initial conditions
    p['traj.phase0.states:x'] = phase0.interp('x', [0,0])
    p['traj.phase0.states:x'] = phase0.interp('x', [0,0])
    p['traj.phase0.states:v'] = phase0.interp('v', [1,-1])
    p['traj.phase0.t_duration'] = 1
    p['traj.phase0.t_initial'] = 0
    
    dm.run_problem(p, make_plots=True)
    

    state history for v state history for x control history for u