From f0daa4f4ac3fd2169d0f2331ac845c2c2332dfde Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Tue, 13 Dec 2022 09:42:39 -0600 Subject: [PATCH 1/7] Added problem from Hull optimal control book --- dymos/examples/hull_problem/__init__.py | 0 dymos/examples/hull_problem/hull_ode.py | 29 ++++++ dymos/examples/hull_problem/test/__init__.py | 0 .../hull_problem/test/test_hull_problem.py | 93 +++++++++++++++++++ .../test/test_hyper_sensitive.py | 2 +- 5 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 dymos/examples/hull_problem/__init__.py create mode 100644 dymos/examples/hull_problem/hull_ode.py create mode 100644 dymos/examples/hull_problem/test/__init__.py create mode 100644 dymos/examples/hull_problem/test/test_hull_problem.py diff --git a/dymos/examples/hull_problem/__init__.py b/dymos/examples/hull_problem/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/dymos/examples/hull_problem/hull_ode.py b/dymos/examples/hull_problem/hull_ode.py new file mode 100644 index 000000000..e7ce08b8a --- /dev/null +++ b/dymos/examples/hull_problem/hull_ode.py @@ -0,0 +1,29 @@ +import numpy as np +import openmdao.api as om + + +class HullProblemODE(om.ExplicitComponent): + def initialize(self): + self.options.declare('num_nodes', types=int) + + def setup(self): + nn = self.options['num_nodes'] + + self.add_input('u', val=np.zeros(nn), desc='control') + + self.add_output('x_dot', val=np.zeros(nn), desc='state rate', units='1/s') + self.add_output('L', val=np.zeros(nn), desc='Lagrangian', units='1/s') + + # Setup partials + self.declare_partials(of='x_dot', wrt='u', rows=np.arange(nn), cols=np.arange(nn), val=1) + self.declare_partials(of='L', wrt='u', rows=np.arange(nn), cols=np.arange(nn)) + + def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): + u = inputs['u'] + + outputs['x_dot'] = u + outputs['L'] = 0.5 * u ** 2 + + def compute_partials(self, inputs, partials, discrete_inputs=None): + u = inputs['u'] + partials['L', 'u'] = u diff --git a/dymos/examples/hull_problem/test/__init__.py b/dymos/examples/hull_problem/test/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/dymos/examples/hull_problem/test/test_hull_problem.py b/dymos/examples/hull_problem/test/test_hull_problem.py new file mode 100644 index 000000000..678c8934b --- /dev/null +++ b/dymos/examples/hull_problem/test/test_hull_problem.py @@ -0,0 +1,93 @@ +from __future__ import print_function, division, absolute_import + +import os +import unittest +import warnings + +import numpy as np + +from openmdao.api import Problem, Group, pyOptSparseDriver, ExecComp +from openmdao.utils.assert_utils import assert_near_equal +from openmdao.utils.general_utils import printoptions +from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse + +import dymos as dm +from dymos import Trajectory, GaussLobatto, Phase, Radau, ExplicitShooting +from dymos.examples.hull_problem.hull_ode import HullProblemODE + +c = 5 + + +@use_tempdirs +class TestHyperSensitive(unittest.TestCase): + + @staticmethod + def make_problem(transcription=GaussLobatto, numseg=30): + p = Problem(model=Group()) + p.driver = pyOptSparseDriver() + p.driver.declare_coloring() + p.driver.options['optimizer'] = 'IPOPT' + p.driver.opt_settings['hessian_approximation'] = 'limited-memory' + p.driver.opt_settings['print_level'] = 0 + p.driver.opt_settings['linear_solver'] = 'mumps' + p.driver.declare_coloring() + + traj = p.model.add_subsystem('traj', Trajectory()) + phase0 = traj.add_phase('phase', Phase(ode_class=HullProblemODE, + transcription=transcription(num_segments=numseg, order=3))) + phase0.set_time_options(fix_initial=True, fix_duration=True) + phase0.add_state('x', fix_initial=True, fix_final=False, rate_source='x_dot') + phase0.add_state('xL', fix_initial=True, fix_final=False, rate_source='L') + phase0.add_control('u', opt=True, targets=['u']) + + obj = ExecComp(f'obj = {c}*x**2/2 + xL') + p.model.add_subsystem('obj_comp', obj) + p.model.add_objective('obj_comp.obj') + p.model.connect('traj.phase.states:x', 'obj_comp.x', src_indices=[-1]) + p.model.connect('traj.phase.states:xL', 'obj_comp.xL', src_indices=[-1]) + + p.setup(check=True) + + p.set_val('traj.phase.states:x', phase0.interp('x', [1.5, 1])) + p.set_val('traj.phase.states:xL', phase0.interp('xL', [0, 1])) + p.set_val('traj.phase.t_initial', 0) + p.set_val('traj.phase.t_duration', 10) + p.set_val('traj.phase.controls:u', phase0.interp('u', [-7, -0.14])) + + return p + + @staticmethod + def solution(x0, td): + + xf = x0 - c * x0 * td / (1 + c * td) + uf = -c * x0 / (1 + c * td) + + return xf, uf + + def test_hull_gauss_lobatto(self): + p = self.make_problem(transcription=GaussLobatto) + dm.run_problem(p) + + xf, uf = self.solution(1.5, 10) + + assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], + xf, + tolerance=1e-4) + + assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], + uf, + tolerance=1e-4) + + def test_hull_radau(self): + p = self.make_problem(transcription=Radau) + dm.run_problem(p) + + xf, uf = self.solution(1.5, 10) + + assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], + xf, + tolerance=1e-4) + + assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], + uf, + tolerance=1e-4) diff --git a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py index 6c4298bdd..9231298f3 100644 --- a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py +++ b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py @@ -16,7 +16,7 @@ from dymos.examples.hyper_sensitive.hyper_sensitive_ode import HyperSensitiveODE -tf = np.float128(10) +tf = 10 @use_tempdirs From 09d4780c64f711a6f1204ebc913f55365ff0c84b Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Tue, 13 Dec 2022 11:01:19 -0600 Subject: [PATCH 2/7] Modified add_objective to accept expressions --- .../hull_problem/test/test_hull_problem.py | 31 +++++++++++++------ dymos/phase/phase.py | 28 ++++++++++++++--- .../explicit_shooting/explicit_shooting.py | 5 +++ dymos/transcriptions/transcription_base.py | 2 +- 4 files changed, 51 insertions(+), 15 deletions(-) diff --git a/dymos/examples/hull_problem/test/test_hull_problem.py b/dymos/examples/hull_problem/test/test_hull_problem.py index 678c8934b..3effb95e2 100644 --- a/dymos/examples/hull_problem/test/test_hull_problem.py +++ b/dymos/examples/hull_problem/test/test_hull_problem.py @@ -19,9 +19,10 @@ @use_tempdirs -class TestHyperSensitive(unittest.TestCase): +class TestHull(unittest.TestCase): @staticmethod + @require_pyoptsparse(optimizer='IPOPT') def make_problem(transcription=GaussLobatto, numseg=30): p = Problem(model=Group()) p.driver = pyOptSparseDriver() @@ -40,11 +41,7 @@ def make_problem(transcription=GaussLobatto, numseg=30): phase0.add_state('xL', fix_initial=True, fix_final=False, rate_source='L') phase0.add_control('u', opt=True, targets=['u']) - obj = ExecComp(f'obj = {c}*x**2/2 + xL') - p.model.add_subsystem('obj_comp', obj) - p.model.add_objective('obj_comp.obj') - p.model.connect('traj.phase.states:x', 'obj_comp.x', src_indices=[-1]) - p.model.connect('traj.phase.states:xL', 'obj_comp.xL', src_indices=[-1]) + phase0.add_objective(f'J = {c}*x**2/2 + xL') p.setup(check=True) @@ -72,11 +69,11 @@ def test_hull_gauss_lobatto(self): assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], xf, - tolerance=1e-4) + tolerance=1e-6) assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], uf, - tolerance=1e-4) + tolerance=1e-6) def test_hull_radau(self): p = self.make_problem(transcription=Radau) @@ -86,8 +83,22 @@ def test_hull_radau(self): assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], xf, - tolerance=1e-4) + tolerance=1e-6) assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], uf, - tolerance=1e-4) + tolerance=1e-6) + + def test_hull_shooting(self): + p = self.make_problem(transcription=ExplicitShooting) + dm.run_problem(p) + + xf, uf = self.solution(1.5, 10) + + assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], + xf, + tolerance=1e-6) + + assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], + uf, + tolerance=1e-6) diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index b25cd5dd6..a69d28109 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1493,7 +1493,7 @@ def add_timeseries(self, name, transcription, subset='all'): 'subset': subset, 'outputs': {}} - def add_objective(self, name, loc='final', index=None, shape=(1,), ref=None, ref0=None, + def add_objective(self, name, loc='final', index=None, shape=(1,), units=None, ref=None, ref0=None, adder=None, scaler=None, parallel_deriv_color=None): """ Add an objective in the phase. @@ -1514,6 +1514,9 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), ref=None, ref used as the objective, assuming C-ordered flattening. shape : int, optional The shape of the objective variable, at a point in time. + units : str, optional + The units of the objective function. If None, use the units associated with the target. + If provided, must be compatible with the target units. ref : float or ndarray, optional Value of response variable that scales to 1.0 in the driver. ref0 : float or ndarray, optional @@ -1528,15 +1531,32 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), ref=None, ref If specified, this design var will be grouped for parallel derivative calculations with other variables sharing the same parallel_deriv_color. """ - obj_dict = {'loc': loc, + expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@'] + if '=' in name: + is_expr = True + elif '=' not in name and any(opr in name for opr in expr_operators): + raise ValueError(f'The expression provided `{name}` has invalid format. ' + 'Expression may be a single variable or an equation ' + 'of the form `constraint_name = func(vars)`') + else: + is_expr = False + + obj_name = name.split('=')[0].strip() if is_expr else name + + obj_dict = {'name': name, + 'loc': loc, 'index': index, 'shape': shape, + 'units': units, 'ref': ref, 'ref0': ref0, 'adder': adder, 'scaler': scaler, - 'parallel_deriv_color': parallel_deriv_color} - self._objectives[name] = obj_dict + 'parallel_deriv_color': parallel_deriv_color, + 'is_expr': is_expr} + self._objectives[obj_name] = obj_dict + if is_expr and obj_name not in self._timeseries['timeseries']['outputs']: + self.add_timeseries_output(name, output_name=obj_name, units=units, shape=shape) def set_time_options(self, units=_unspecified, fix_initial=_unspecified, fix_duration=_unspecified, input_initial=_unspecified, diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py index 380bec2c8..e334d488c 100644 --- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py +++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py @@ -547,6 +547,11 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None): units = control_rate_units linear = False obj_path = f'polynomial_control_rates:{var}' + elif var_type == 'timeseries_exec_comp_output': + shape = (1,) + units = None + obj_path = f'timeseries.timeseries_exec_comp.{var}' + linear = False else: # Failed to find variable, assume it is in the ODE. This requires introspection. raise NotImplementedError('cannot yet constrain/optimize an ODE output using explicit shooting') diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index 28c7af713..47b90234f 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -535,7 +535,7 @@ def configure_objective(self, phase): index = options['index'] loc = options['loc'] - obj_path, shape, units, _ = self._get_objective_src(name, loc, phase) + obj_path, shape, _, _ = self._get_objective_src(name, loc, phase) shape = options['shape'] if shape is None else shape From 5921823145466d01d32acabe884d351c725f1aae Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Tue, 13 Dec 2022 15:09:26 -0600 Subject: [PATCH 3/7] Added documentation --- .../examples/hull/hull_problem.ipynb | 515 ++++++++++++++++++ .../features/phases/objective.ipynb | 20 +- dymos/examples/hull_problem/hull_ode.py | 3 - .../hull_problem/test/test_hull_problem.py | 10 +- dymos/phase/phase.py | 3 +- 5 files changed, 537 insertions(+), 14 deletions(-) create mode 100644 docs/dymos_book/examples/hull/hull_problem.ipynb diff --git a/docs/dymos_book/examples/hull/hull_problem.ipynb b/docs/dymos_book/examples/hull/hull_problem.ipynb new file mode 100644 index 000000000..d358a1ac6 --- /dev/null +++ b/docs/dymos_book/examples/hull/hull_problem.ipynb @@ -0,0 +1,515 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "e9bef506", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "e9bef506", + "outputId": "23e61fa9-803c-416c-c2ea-49fda74021d7", + "tags": [ + "active-ipynb", + "remove-input", + "remove-output" + ] + }, + "outputs": [], + "source": [ + "# This cell is mandatory in all Dymos documentation notebooks.\n", + "missing_packages = []\n", + "try:\n", + " import openmdao.api as om\n", + "except ImportError:\n", + " if 'google.colab' in str(get_ipython()):\n", + " !python -m pip install openmdao[notebooks]\n", + " else:\n", + " missing_packages.append('openmdao')\n", + "try:\n", + " import dymos as dm\n", + "except ImportError:\n", + " if 'google.colab' in str(get_ipython()):\n", + " !python -m pip install dymos\n", + " else:\n", + " missing_packages.append('dymos')\n", + "try:\n", + " import pyoptsparse\n", + "except ImportError:\n", + " if 'google.colab' in str(get_ipython()):\n", + " !pip install -q condacolab\n", + " import condacolab\n", + " condacolab.install_miniconda()\n", + " !conda install -c conda-forge pyoptsparse\n", + " else:\n", + " missing_packages.append('pyoptsparse')\n", + "if missing_packages:\n", + " raise EnvironmentError('This notebook requires the following packages '\n", + " 'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')" + ] + }, + { + "cell_type": "markdown", + "id": "b0a95015", + "metadata": { + "id": "b0a95015" + }, + "source": [ + "(examples:hull)=\n", + "# The Hull Problem\n", + "\n", + "The Hull problem is a 1-DOF optimal control problem {cite}`hull2003oct`. It can be stated as:\n", + "\n", + "Minimize the control effort required to move a frictionless sliding block from some initial position such that the final displacement from a pre-specified point is minimized." + ] + }, + { + "cell_type": "markdown", + "id": "977155e8", + "metadata": { + "id": "977155e8" + }, + "source": [ + "## State and control variables\n", + "\n", + "This system has one state variables, the position ($x$) of the sliding block. \n", + "\n", + "This system has a single control variable ($u$), the velocity of the block.\n", + "\n", + "The dynamics of the system are governed by\n", + "\n", + "\\begin{align}\n", + " \\dot{x} &= u\n", + "\\end{align}\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "6e236c20", + "metadata": { + "id": "6e236c20" + }, + "source": [ + "## Problem Definition\n", + "\n", + "We seek to minimize the control effort required and minimize the displacement from the origin.\n", + "\n", + "\\begin{align}\n", + " \\mathrm{Minimize} \\, J &= 2.5x_f^2 \\, + \\, 0.5 \\int_0^1 u^2 dt\n", + "\\end{align}\n", + "\n", + "Subject to the initial conditions\n", + "\n", + "\\begin{align}\n", + " t_0 &= 0.0 \\\\\n", + " x_0 &= 1.5\n", + "\\end{align}\n", + "\n", + "and the terminal constraints\n", + "\n", + "\\begin{align}\n", + " t_f &= 10.0\n", + "\\end{align}\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "1229fca5", + "metadata": {}, + "source": [ + "## Dealing with combined terminal and integral costs in Dymos\n", + "\n", + "In classic optimal control, the objective is often broken into the terminal component (the Mayer term) and the integral component (the Lagrange term).\n", + "Dymos does not distinguish between the two.\n", + "In this case, since the objective $J$ consists of both a terminal cost and an integrated cost (Bolza form), we add a term to the ODE to account for the integrated quantity\n", + "\n", + "\n", + "\\begin{align}\n", + " \\dot{x_L} &= L \\\\\n", + " L &= 0.5 u^2\n", + "\\end{align}\n", + "\n", + "where $x_L$ is a state added to account for the Lagrange term.\n", + "\n", + "Dymos supports the definition of simple mathematical expressions as the cost, so the final value of $x_L$ can be added to the final value of $2.5x^2$." + ] + }, + { + "cell_type": "markdown", + "id": "03ef4a09", + "metadata": { + "id": "03ef4a09" + }, + "source": [ + "## Defining the ODE\n", + "\n", + "The following code implements the equations of motion for the Hull problem.\n", + "Since the rate of $x$ is given by a control ($u$), there is no need to compute its rate in the ODE.\n", + "Dymos can pull their values from those other states and controls.\n", + "The ODE, therefore, only needs to compute the rate of change of $x_L$ ($L$).\n", + "\n", + "A few things to note:\n", + "\n", + "1. By providing the tag `dymos.state_rate_source:{name}`, we're letting Dymos know what states need to be integrated, there's no need to specify a rate source when using this ODE in our Phase.\n", + "2. Pairing the above tag with `dymos.state_units:{units}` means we don't have to specify units when setting properties for the state in our run script." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec2b8a83", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "ec2b8a83", + "outputId": "668bf949-c304-4e34-9fca-27697a497d4b" + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import openmdao.api as om\n", + "\n", + "\n", + "class HullProblemODE(om.ExplicitComponent):\n", + " def initialize(self):\n", + " self.options.declare('num_nodes', types=int)\n", + "\n", + " def setup(self):\n", + " nn = self.options['num_nodes']\n", + "\n", + " self.add_input('u', val=np.zeros(nn), desc='control')\n", + "\n", + " self.add_output('L', val=np.zeros(nn), desc='Lagrangian', units='1/s')\n", + "\n", + " # Setup partials\n", + " self.declare_partials(of='L', wrt='u', rows=np.arange(nn), cols=np.arange(nn))\n", + "\n", + " def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):\n", + " u = inputs['u']\n", + "\n", + " outputs['L'] = 0.5 * u ** 2\n", + "\n", + " def compute_partials(self, inputs, partials, discrete_inputs=None):\n", + " u = inputs['u']\n", + " partials['L', 'u'] = u\n" + ] + }, + { + "cell_type": "markdown", + "id": "0031d62b", + "metadata": { + "id": "hVK50KxH6YJ4" + }, + "source": [ + "## Solving the Hull problem with Dymos\n", + "\n", + "The following script solves the Hull problem with Dymos.\n", + "\n", + "To begin, import the packages we require:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ii-ApZna669K", + "metadata": { + "id": "ii-ApZna669K" + }, + "outputs": [], + "source": [ + "import dymos as dm\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "OI5v_AkL7qYm", + "metadata": { + "id": "OI5v_AkL7qYm" + }, + "source": [ + "We then instantiate an OpenMDAO problem and set the optimizer and its options.\n", + "\n", + "The call to `declare_coloring` tells the optimizer to attempt to find a sparsity pattern that minimizes the work required to compute the derivatives across the model.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3yXzgcIw8Tjq", + "metadata": { + "id": "3yXzgcIw8Tjq" + }, + "outputs": [], + "source": [ + "#\n", + "# Initialize the Problem and the optimization driver\n", + "#\n", + "p = om.Problem()\n", + " \n", + "p.driver = om.pyOptSparseDriver()\n", + "p.driver.declare_coloring()" + ] + }, + { + "cell_type": "markdown", + "id": "OcugnHOL8fIF", + "metadata": { + "id": "OcugnHOL8fIF" + }, + "source": [ + "Next, we add a Dymos Trajectory group to the problem's model and add a phase to it.\n", + "\n", + "In this case we're using the Radau pseudospectral transcription to solve the problem." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "sGYfpaiI8eH-", + "metadata": { + "id": "sGYfpaiI8eH-" + }, + "outputs": [], + "source": [ + "#\n", + "# Create a trajectory and add a phase to it\n", + "#\n", + "traj = p.model.add_subsystem('traj', dm.Trajectory())\n", + "tx = transcription=dm.Radau(num_segments=24)\n", + "phase = traj.add_phase('phase0', dm.Phase(ode_class=HullProblemODE, transcription=tx))" + ] + }, + { + "cell_type": "markdown", + "id": "IjCxNJQV82u8", + "metadata": { + "id": "IjCxNJQV82u8" + }, + "source": [ + "At this point, we set the options on the main variables used in a Dymos phase. \n", + "\n", + "In addition to `time`, we have two states (`x`, and `x_L`) and a single control (`u`). \n", + "\n", + "Here we use bounds on the states themselves to constrain the initial values of `x` and1 `x_L`.\n", + "From an optimization perspective, this means that we are removing the first and last values in the state histories of $x$ and $x_L$ from the vector of design variables.\n", + "Their initial values will remain unchanged throughout the optimization process.\n", + "\n", + "On the other hand, we could specify `fix_initial=False, fix_final=False` for these values, and Dymos would be free to change them.\n", + "We would then need to put a boundary constraint in place to enforce their final values.\n", + "Feel free to experiment with different ways of enforcing the boundary constraints on this problem and see how it affects performance.\n", + "\n", + "The scaler values (`ref`) are all set to 1 here.\n", + "\n", + "Also, we don't need to specify targets for any of the variables here because their names _are_ the targets in the top-level of the model.\n", + "The rate source and units for the states are obtained from the tags in the ODE component we previously defined.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "IdjlON_895PU", + "metadata": { + "id": "IdjlON_895PU" + }, + "outputs": [], + "source": [ + "#\n", + "# Set the variables\n", + "#\n", + "phase.set_time_options(fix_initial=True, fix_duration=True)\n", + "\n", + "phase.set_time_options(fix_initial=True, fix_duration=True)\n", + "phase.add_state('x', fix_initial=True, fix_final=False, rate_source='u')\n", + "phase.add_state('xL', fix_initial=True, fix_final=False, rate_source='L')\n", + "phase.add_control('u', opt=True, targets=['u'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d498cb7b", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "d498cb7b", + "outputId": "18be63ac-a3cb-4632-eb7d-130eced98d70" + }, + "outputs": [], + "source": [ + "#\n", + "# Minimize time at the end of the phase\n", + "#\n", + "phase.add_objective('J = 2.5*x**2 + xL')\n", + "\n", + "#\n", + "# Setup the Problem\n", + "#\n", + "p.setup()" + ] + }, + { + "cell_type": "markdown", + "id": "E6q4fW_NAx6B", + "metadata": { + "id": "E6q4fW_NAx6B" + }, + "source": [ + "We then set the initial guesses for the variables in the problem and solve it.\n", + "\n", + "We're using the phase `interp` method to provide initial guesses for the states and controls.\n", + "In this case, by giving it two values, it is linearly interpolating from the first value to the second value, and then returning the interpolated value at the input nodes for the given variable.\n", + "\n", + "Finally, we use the `dymos.run_problem` method to execute the problem.\n", + "This interface allows us to do some things that the standard OpenMDAO `problem.run_driver` interface does not.\n", + "It will automatically record the final solution achieved by the optimizer in case named `'final'` in a file called `dymos_solution.db`.\n", + "By specifying `simulate=True`, it will automatically follow the solution with an explicit integration using `scipy.solve_ivp`.\n", + "The results of the simulation are stored in a case named `final` in the file `dymos_simulation.db`.\n", + "This explicit simulation demonstrates how the system evolved with the given controls, and serves as a check that we're using a dense enough grid (enough segments and segments of sufficient order) to accurately represent the solution.\n", + "\n", + "If those two solution didn't agree reasonably well, we could rerun the problem with a more dense grid.\n", + "Instead, we're asking Dymos to automatically change the grid if necessary by specifying `refine_method='ph'`.\n", + "This will attempt to repeatedly solve the problem and change the number of segments and segment orders until the solution is in reasonable agreement." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "iRY53Rq0_0c6", + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "iRY53Rq0_0c6", + "outputId": "1d987fbe-e703-4e89-cb8c-14af763b8d7f", + "tags": [ + "hide-output" + ] + }, + "outputs": [], + "source": [ + "#\n", + "# Set the initial values\n", + "#\n", + "p.set_val('traj.phase0.states:x', phase.interp('x', [1.5, 1]))\n", + "p.set_val('traj.phase0.states:xL', phase.interp('xL', [0, 1]))\n", + "p.set_val('traj.phase0.t_initial', 0)\n", + "p.set_val('traj.phase0.t_duration', 10)\n", + "p.set_val('traj.phase0.controls:u', phase.interp('u', [-7, -0.14]))\n", + "#\n", + "# Solve for the optimal trajectory\n", + "#\n", + "dm.run_problem(p, run_driver=True, simulate=True)\n" + ] + }, + { + "cell_type": "markdown", + "id": "3f9a5d44", + "metadata": { + "id": "3f9a5d44" + }, + "source": [ + "## Plotting the solution\n", + "\n", + "The recommended practice is to obtain values from the recorded cases.\n", + "While the problem object can also be queried for values, building plotting scripts that use the case recorder files as the data source means that the problem doesn't need to be solved just to change a plot.\n", + "Here we load values of various variables from the solution and simulation for use in the animation to follow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "beea4c6a", + "metadata": { + "id": "beea4c6a", + "scrolled": true + }, + "outputs": [], + "source": [ + "sol = om.CaseReader('dymos_solution.db').get_case('final')\n", + "sim = om.CaseReader('dymos_simulation.db').get_case('final')\n", + "\n", + "t = sol.get_val('traj.phase0.timeseries.time')\n", + "x = sol.get_val('traj.phase0.timeseries.states:x')\n", + "xL = sol.get_val('traj.phase0.timeseries.states:xL')\n", + "u = sol.get_val('traj.phase0.timeseries.controls:u')\n", + "\n", + "t_sim = sim.get_val('traj.phase0.timeseries.time')\n", + "x_sim = sim.get_val('traj.phase0.timeseries.states:x')\n", + "xL_sim = sim.get_val('traj.phase0.timeseries.states:xL')\n", + "u_sim = sim.get_val('traj.phase0.timeseries.controls:u')\n", + "\n", + "fig = plt.figure(constrained_layout=True, figsize=(12, 4))\n", + "gs = fig.add_gridspec(3, 1)\n", + "\n", + "x_ax = fig.add_subplot(gs[0, 0])\n", + "xL_ax = fig.add_subplot(gs[1, 0])\n", + "u_ax = fig.add_subplot(gs[2, 0])\n", + "\n", + "x_ax.set_ylabel('x ($m$)')\n", + "xL_ax.set_ylabel('xL ($m/s$)')\n", + "u_ax.set_ylabel('u ($m/s^2$)')\n", + "xL_ax.set_xlabel('t (s)')\n", + "u_ax.set_xlabel('t (s)')\n", + "\n", + "x_sol_handle, = x_ax.plot(t, x, 'o', ms=1)\n", + "xL_ax.plot(t, xL, 'o', ms=1)\n", + "u_ax.plot(t, u, 'o', ms=1)\n", + "\n", + "x_sim_handle, = x_ax.plot(t_sim, x_sim, '-', ms=1)\n", + "xL_ax.plot(t_sim, xL_sim, '-', ms=1)\n", + "u_ax.plot(t_sim, u_sim, '-', ms=1)\n", + "\n", + "for ax in [x_ax, xL_ax, u_ax]:\n", + " ax.grid(True, alpha=0.2)\n", + " \n", + "plt.figlegend([x_sol_handle, x_sim_handle], ['solution', 'simulation'], ncol=2, loc='lower center');\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "6d27dc72", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "```{bibliography}\n", + ":filter: docname in docnames\n", + "```" + ] + } + ], + "metadata": { + "celltoolbar": "Tags", + "colab": { + "collapsed_sections": [], + "name": "mountain_car.ipynb", + "provenance": [] + }, + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/dymos_book/features/phases/objective.ipynb b/docs/dymos_book/features/phases/objective.ipynb index c6a39a636..82dcabdc0 100644 --- a/docs/dymos_book/features/phases/objective.ipynb +++ b/docs/dymos_book/features/phases/objective.ipynb @@ -66,6 +66,7 @@ "- the name of a control variable\n", "- the name of a control rate or rate2 (second derivative)\n", "- the path of an output in the ODE relative to the top level of the ODE\n", + "- an expression that is a combination of any of the above in the form of an equation\n", "\n", "Dymos will find the full path of the given variable and add it to the problem as an objective.\n", "\n", @@ -95,8 +96,23 @@ "\n", "```python\n", "phase.add_objective('mass', loc='final', scaler=-1)\n", + "```\n", + "\n", + "## Example: Minimizing Final Displacement\n", + "\n", + "This example assumes that the phase has a state variable named *x*. Here, *disp* is added to the timeseries and the final value computed is minimized\n", + "\n", + "```python\n", + "phase.add_objective('disp=x**2', loc='final')\n", "```" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -109,7 +125,7 @@ } }, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -123,7 +139,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.1" + "version": "3.9.12" } }, "nbformat": 4, diff --git a/dymos/examples/hull_problem/hull_ode.py b/dymos/examples/hull_problem/hull_ode.py index e7ce08b8a..0fbdef997 100644 --- a/dymos/examples/hull_problem/hull_ode.py +++ b/dymos/examples/hull_problem/hull_ode.py @@ -11,17 +11,14 @@ def setup(self): self.add_input('u', val=np.zeros(nn), desc='control') - self.add_output('x_dot', val=np.zeros(nn), desc='state rate', units='1/s') self.add_output('L', val=np.zeros(nn), desc='Lagrangian', units='1/s') # Setup partials - self.declare_partials(of='x_dot', wrt='u', rows=np.arange(nn), cols=np.arange(nn), val=1) self.declare_partials(of='L', wrt='u', rows=np.arange(nn), cols=np.arange(nn)) def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): u = inputs['u'] - outputs['x_dot'] = u outputs['L'] = 0.5 * u ** 2 def compute_partials(self, inputs, partials, discrete_inputs=None): diff --git a/dymos/examples/hull_problem/test/test_hull_problem.py b/dymos/examples/hull_problem/test/test_hull_problem.py index 3effb95e2..bc040f5a6 100644 --- a/dymos/examples/hull_problem/test/test_hull_problem.py +++ b/dymos/examples/hull_problem/test/test_hull_problem.py @@ -1,14 +1,8 @@ from __future__ import print_function, division, absolute_import - -import os import unittest -import warnings - -import numpy as np -from openmdao.api import Problem, Group, pyOptSparseDriver, ExecComp +from openmdao.api import Problem, Group, pyOptSparseDriver from openmdao.utils.assert_utils import assert_near_equal -from openmdao.utils.general_utils import printoptions from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse import dymos as dm @@ -37,7 +31,7 @@ def make_problem(transcription=GaussLobatto, numseg=30): phase0 = traj.add_phase('phase', Phase(ode_class=HullProblemODE, transcription=transcription(num_segments=numseg, order=3))) phase0.set_time_options(fix_initial=True, fix_duration=True) - phase0.add_state('x', fix_initial=True, fix_final=False, rate_source='x_dot') + phase0.add_state('x', fix_initial=True, fix_final=False, rate_source='u') phase0.add_state('xL', fix_initial=True, fix_final=False, rate_source='L') phase0.add_control('u', opt=True, targets=['u']) diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index a69d28109..b4e386dc5 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -1505,7 +1505,8 @@ def add_objective(self, name, loc='final', index=None, shape=(1,), units=None, r ---------- name : str Name of the objective variable. This should be one of the integration variable, a state or control - variable, or the path to an output from the top level of the RHS. + variable, the path to an output from the top level of the RHS, or an expression to be evaluated. + If an expression, it must be provided in the form of an equation with a left- and right-hand side. loc : str Where in the phase the objective is to be evaluated. Valid options are 'initial' and 'final'. The default is 'final'. From 8669aa9529359eb57d2a4b4179acff5a8e24d9d3 Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Tue, 13 Dec 2022 15:13:53 -0600 Subject: [PATCH 4/7] Updated bibliography --- docs/dymos_book/references.bib | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/dymos_book/references.bib b/docs/dymos_book/references.bib index bd2fac251..381796caa 100644 --- a/docs/dymos_book/references.bib +++ b/docs/dymos_book/references.bib @@ -588,3 +588,10 @@ @article{Kelly2017 issn = {00361445}, keywords = {Direct collocation, Direct transcription, Optimal control, Robotics, Trajectory optimization, Tutorial} } + +@book{hull2003oct, + title={Optimal Control Theory for Applications}, + author={Hull, David G}, + year={2003}, + publisher={Springer} +} From 9e306320ff93f925afbb767bb4053417df9cdf50 Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Tue, 13 Dec 2022 15:41:51 -0600 Subject: [PATCH 5/7] Added hull file to toc --- docs/dymos_book/_toc.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml index 638b90387..29561f73f 100644 --- a/docs/dymos_book/_toc.yml +++ b/docs/dymos_book/_toc.yml @@ -24,6 +24,7 @@ parts: - file: examples/bryson_denham/bryson_denham - file: examples/commercial_aircraft/commercial_aircraft - file: examples/double_integrator/double_integrator + - file: examples/hull/hull_problem - file: examples/hypersensitive/hypersensitive - file: examples/finite_burn_orbit_raise/finite_burn_orbit_raise - file: examples/length_constrained_brachistochrone/length_constrained_brachistochrone From 35da2844ea84692e744df2b75c07c81f61b0b5be Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Wed, 14 Dec 2022 07:38:05 -0600 Subject: [PATCH 6/7] Fixed a bug in ODE evaluation group for shooting --- .../transcriptions/explicit_shooting/ode_evaluation_group.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py index cc098e797..c452d4769 100644 --- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py +++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py @@ -381,10 +381,10 @@ def _get_rate_source_path(self, state_var): rate_path = f'states:{var}' io = 'input' elif self.control_options is not None and var in self.control_options: - rate_path = f'controls:{var}' + rate_path = f'control_values:{var}' io = 'output' elif self.polynomial_control_options is not None and var in self.polynomial_control_options: - rate_path = f'polynomial_controls:{var}' + rate_path = f'polynomial_control_values:{var}' io = 'output' elif self.parameter_options is not None and var in self.parameter_options: rate_path = f'parameter_vals:{var}' From cfe4050a236bfac63b3c69eabc4d2890b579e084 Mon Sep 17 00:00:00 2001 From: Kaushik Ponnapalli Date: Wed, 14 Dec 2022 16:50:18 -0600 Subject: [PATCH 7/7] loosened tolerances on hull problem --- .../examples/hull_problem/test/test_hull_problem.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/dymos/examples/hull_problem/test/test_hull_problem.py b/dymos/examples/hull_problem/test/test_hull_problem.py index bc040f5a6..a9734058c 100644 --- a/dymos/examples/hull_problem/test/test_hull_problem.py +++ b/dymos/examples/hull_problem/test/test_hull_problem.py @@ -63,11 +63,11 @@ def test_hull_gauss_lobatto(self): assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], xf, - tolerance=1e-6) + tolerance=1e-4) assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], uf, - tolerance=1e-6) + tolerance=1e-4) def test_hull_radau(self): p = self.make_problem(transcription=Radau) @@ -77,11 +77,11 @@ def test_hull_radau(self): assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], xf, - tolerance=1e-6) + tolerance=1e-4) assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], uf, - tolerance=1e-6) + tolerance=1e-4) def test_hull_shooting(self): p = self.make_problem(transcription=ExplicitShooting) @@ -91,8 +91,8 @@ def test_hull_shooting(self): assert_near_equal(p.get_val('traj.phase.timeseries.states:x')[-1], xf, - tolerance=1e-6) + tolerance=1e-4) assert_near_equal(p.get_val('traj.phase.timeseries.controls:u')[-1], uf, - tolerance=1e-6) + tolerance=1e-4)