.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples\plot_betts_10_144_145_opty.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_plot_betts_10_144_145_opty.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_plot_betts_10_144_145_opty.py:


Van der Pol Oscillator  (Betts 10.144 / 145)
============================================

https://en.wikipedia.org/wiki/Van_der_Pol_oscillator

These are examples 10.144 / 145 from *John T. Betts, Practical Methods for
Optimal Control Using NonlinearProgramming*, 3rd edition, Chapter 10: Test
Problems.
It is described in more detail in section 4.14. example 4.11 of the book.

Note
----

- A rather high number of nodes is required to get a good solution. Maybe
  an indication that the problem is stiff.
- Both formulations give very similar results.

.. GENERATED FROM PYTHON SOURCE LINES 21-29

.. code-block:: Python

    import numpy as np
    import sympy as sm
    import sympy.physics.mechanics as me
    import time
    from opty.direct_collocation import Problem
    from opty.utils import create_objective_function, MathJaxRepr
    import matplotlib.pyplot as plt








.. GENERATED FROM PYTHON SOURCE LINES 30-43

First version of the problem (10.144)
=====================================

**States**

:math:`y_1, y_2` : state variables

**Controls**

:math:`u` : control variables

Equations of Motion.
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 43-55

.. code-block:: Python

    t = me.dynamicsymbols._t

    y1, y2 = me.dynamicsymbols('y1, y2')
    u = me.dynamicsymbols('u')

    eom = sm.Matrix([
        -y1.diff() + y2,
        -y2.diff(t) + (1 - y1**2)*y2 - y1 + u,
    ])

    MathJaxRepr(eom)






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    $$\left[\begin{matrix}y_{2} - \dot{y}_{1}\\\left(1 - y_{1}^{2}\right) y_{2} + u - y_{1} - \dot{y}_{2}\end{matrix}\right]$$
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 56-58

Define and Solve the Optimization Problem.
------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 58-98

.. code-block:: Python

    num_nodes = 20001

    t0, tf = 0.0, 5.0
    interval_value = (tf - t0) / (num_nodes - 1)

    state_symbols = (y1, y2)
    unkonwn_input_trajectories = (u, )

    objective = sm.integrate(u**2 + y1**2 + y2**2, t)
    obj, obj_grad = create_objective_function(
        objective,
        state_symbols,
        unkonwn_input_trajectories,
        tuple(),
        num_nodes,
        interval_value,
        time_symbol=t
    )

    instance_constraints = (
            y1.func(t0) - 1.0,
            y2.func(t0),
    )

    bounds = {
            y2: (-0.4, np.inf),
    }

    prob = Problem(
        obj,
        obj_grad,
        eom,
        state_symbols,
        num_nodes,
        interval_value,
        instance_constraints=instance_constraints,
        bounds=bounds,
        time_symbol=t,
    )








.. GENERATED FROM PYTHON SOURCE LINES 99-101

Solve the optimization problem.
Give some rough estimates for the trajectories.

.. GENERATED FROM PYTHON SOURCE LINES 101-104

.. code-block:: Python


    initial_guess = np.ones(prob.num_free)








.. GENERATED FROM PYTHON SOURCE LINES 105-106

Find the optimal solution.

.. GENERATED FROM PYTHON SOURCE LINES 106-116

.. code-block:: Python

    start = time.time()
    solution, info = prob.solve(initial_guess)
    end = time.time()
    print(f"Solving took {end - start:.2f} seconds.")
    print(info['status_msg'])
    Jstar = 2.95369916
    print(f"Objectve is: {info['obj_val']:.8f}, " +
          f"as per the book it is {Jstar}, so the deviation is: "
          f"{(info['obj_val'] - Jstar) / Jstar * 100:.5e} %")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Solving took 4.67 seconds.
    b'Algorithm terminated successfully at a locally optimal point, satisfying the convergence tolerances (can be specified by options).'
    Objectve is: 2.95319329, as per the book it is 2.95369916, so the deviation is: -1.71266e-02 %




.. GENERATED FROM PYTHON SOURCE LINES 117-118

Store the first solution for later comparison.

.. GENERATED FROM PYTHON SOURCE LINES 118-119

.. code-block:: Python

    solution1 = solution







.. GENERATED FROM PYTHON SOURCE LINES 120-121

Plot the optimal state and input trajectories.

.. GENERATED FROM PYTHON SOURCE LINES 121-123

.. code-block:: Python

    _ = prob.plot_trajectories(solution, show_bounds=True)




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_001.png
   :alt: State Trajectories, Input Trajectories
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 124-125

Plot the constraint violations.

.. GENERATED FROM PYTHON SOURCE LINES 125-127

.. code-block:: Python

    _ = prob.plot_constraint_violations(solution)




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_002.png
   :alt: Constraint violations
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 128-129

Plot the objective function as a function of optimizer iteration.

.. GENERATED FROM PYTHON SOURCE LINES 129-131

.. code-block:: Python

    _ = prob.plot_objective_value()




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_003.png
   :alt: Objective Value
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 132-148

Second version of the problem (10.145)
======================================

It is same as problem 10.144 but reformulated.
It has two control variables and one additional algebraic equation of motion.

**States**

:math:`y_1, y_2` : state variables

**Controls**

:math:`u, v` : control variables

Equations of Motion.
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 148-158

.. code-block:: Python

    y1, y2, v = me.dynamicsymbols('y1, y2, v')
    u = me.dynamicsymbols('u')

    eom = sm.Matrix([
        -y1.diff() + y2,
        -y2.diff(t) + v - y1 + u,
        v - (1-y1**2)*y2,
    ])
    MathJaxRepr(eom)






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    $$\left[\begin{matrix}y_{2} - \dot{y}_{1}\\u + v - y_{1} - \dot{y}_{2}\\- \left(1 - y_{1}^{2}\right) y_{2} + v\end{matrix}\right]$$
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 159-161

Define and Solve the Optimization Problem
-----------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 161-195

.. code-block:: Python

    state_symbols = (y1, y2)
    unkonwn_input_trajectories = (u, v)

    objective = sm.integrate(u**2 + y1**2 + y2**2, t)
    obj, obj_grad = create_objective_function(
        objective,
        state_symbols,
        unkonwn_input_trajectories,
        tuple(),
        num_nodes,
        interval_value
    )

    instance_constraints = (
            y1.func(t0) - 1.0,
            y2.func(t0),
    )

    bounds = {
            y2: (-0.4, np.inf),
    }

    prob = Problem(
        obj,
        obj_grad,
        eom,
        state_symbols,
        num_nodes,
        interval_value,
        instance_constraints=instance_constraints,
        bounds=bounds,
        time_symbol=t,
    )








.. GENERATED FROM PYTHON SOURCE LINES 196-198

Solve the optimization problem.
Give some rough estimates for the trajectories.

.. GENERATED FROM PYTHON SOURCE LINES 198-201

.. code-block:: Python


    initial_guess = np.ones(prob.num_free)








.. GENERATED FROM PYTHON SOURCE LINES 202-203

Find the optimal solution.

.. GENERATED FROM PYTHON SOURCE LINES 203-212

.. code-block:: Python

    start = time.time()
    solution, info = prob.solve(initial_guess)
    end = time.time()
    print(f"Solving took {end - start:.2f} seconds.")
    Jstar = 2.95369916
    print(f"Objectve is: {info['obj_val']:.8f}, " +
          f"as per the book it is {Jstar}, so the deviation is: "
          f"{(info['obj_val'] - Jstar) / Jstar * 100:.5e} %")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Solving took 6.36 seconds.
    Objectve is: 2.95319329, as per the book it is 2.95369916, so the deviation is: -1.71266e-02 %




.. GENERATED FROM PYTHON SOURCE LINES 213-214

Plot the optimal state and input trajectories.

.. GENERATED FROM PYTHON SOURCE LINES 214-216

.. code-block:: Python

    _ = prob.plot_trajectories(solution, show_bounds=True)




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_004.png
   :alt: State Trajectories, Input Trajectories
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 217-218

Plot the constraint violations.

.. GENERATED FROM PYTHON SOURCE LINES 218-220

.. code-block:: Python

    _ = prob.plot_constraint_violations(solution)




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_005.png
   :alt: Constraint violations
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 221-222

Plot the objective function as a function of optimizer iteration.

.. GENERATED FROM PYTHON SOURCE LINES 222-224

.. code-block:: Python

    _ = prob.plot_objective_value()




.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_006.png
   :alt: Objective Value
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 225-227

Plot the Difference between the two Solutions.
----------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 227-242

.. code-block:: Python

    diffy1 = solution1[: num_nodes] - solution[: num_nodes]
    diffy2 = solution1[num_nodes: 2*num_nodes] - solution[num_nodes: 2*num_nodes]
    diffu = solution1[2*num_nodes:] - solution[2*num_nodes: 3*num_nodes]
    times = np.linspace(t0, tf, num_nodes)

    fig, ax = plt.subplots(3, 1, figsize=(7, 4), sharex=True,
                           constrained_layout=True)
    ax[0].plot(times, diffy1, label='Delta y1')
    ax[0].legend()
    ax[1].plot(times, diffy2, label='Delta y2')
    ax[1].legend()
    ax[2].plot(times, diffu, label='Delta u')
    ax[2].legend()
    ax[2].set_xlabel('Time')
    _ = ax[0].set_title('Differences between the two solutions')



.. image-sg:: /examples/images/sphx_glr_plot_betts_10_144_145_opty_007.png
   :alt: Differences between the two solutions
   :srcset: /examples/images/sphx_glr_plot_betts_10_144_145_opty_007.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 27.031 seconds)


.. _sphx_glr_download_examples_plot_betts_10_144_145_opty.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_betts_10_144_145_opty.ipynb <plot_betts_10_144_145_opty.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_betts_10_144_145_opty.py <plot_betts_10_144_145_opty.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_betts_10_144_145_opty.zip <plot_betts_10_144_145_opty.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_