opty does not accept only one scalar differential equation

[Peter230655]

Problem 10.7 in Chapter 10 of John T. Betts' Practical Methods for....., 3rd edition, only has one scalar differential equation: $\dot y = - y^3 + u$, where u is the control. When I run it, I get this error:

When I 'artificially' enlarge the system: $uy = \dfrac{d}{dt} y$ $uy = - y^3 + u$ all works fine and it find the solution. This is definitely not a big thing, it is probably a rare case that there is only one scalar DE, but I wanted to point it out.

[moorepants]

Yes, opty assumes you are working with a physics system that has a mass and forces, thus you always get at least one second order differential equation.

[Peter230655]

Yes, opty assumes you are working with a physics system that has a mass and forces, thus you always get at least one second order differential equation.

Makes sense. Betts gives the dynamic equations as a set of first order DEs, at least in the three examples I have calculated. In example 10.6 it is a set of 10 first order DEs, which do not look like mechanical equations - yet the results look fine. So, as long as it is more than one first order DEs in the eoms, opty does not care how they 'originated?

[moorepants]

So, as long as it is more than one first order DEs in the eoms, opty does not care how they 'originated?

I think that the only assumption is that there are n states and n equations, where n is > 1.

Two things:

• An example in the gallery that shows how to use a single differential equation with your hack is useful as a temporary solution.
• We should relax the assumption opty makes, because we also would like to add extra constraint equations where n states != num constraint equations (e.g. in https://github.com/csu-hmc/opty/issues/156).

[Peter230655]

1. I will push a 'beautified' version of example 10.7 to opty-gallery (it is in pst_notebooks, but combined with 10.8 - an interesting comparison in my opinion in its own right)
2. Should I raise an issue, or is this implicit in #156?

[moorepants]

This issues is sufficient for the request for support for single first order differential equations.

[Peter230655]

This issues is sufficient for the request for support for single first order differential equations.

Not sure I understand. Raise a new issue or not? (Sorry, my command of the English language...)

[moorepants]

This is the issue, no new issue is needed.

[moorepants]

You could try removing https://github.com/csu-hmc/opty/blob/master/opty/direct_collocation.py#L1184 and just see what happens.

[Peter230655]

You could try removing https://github.com/csu-hmc/opty/blob/master/opty/direct_collocation.py#L1184 and just see what happens.

I will try in the next couple of days and report.

[Peter230655]

This program runs fine im VSC, but in examples-gallery it looks like it never evewn starts to be exceuted. Did anything change in examples-gallery lately, or is something wrong with my program below for examples-gallery? Thanks!

# %%
"""
Hypersensitive Control
======================

This is example 10.7 from Betts' Test Problems.

**States**

- y : state variable
- uy: its speed

**Specifieds**

- u : control variable

Note: the state variable uy is needed because opt currently needs minimum two
differential equations in the equations of motion. Mathematically it is not
needed.

"""

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

# %%
# Equations of motion.

t = me.dynamicsymbols._t
y, uy, u = me.dynamicsymbols('y uy, u')

eom = sm.Matrix([ uy - y.diff(t), -uy - y**3 + u])
sm.pprint(eom)

# %%

def solve_optimization(nodes, tf):
    t0, tf = 0.0, tf
    num_nodes = nodes
    interval_value = (tf - t0)/(num_nodes - 1)

    # Provide some reasonably realistic values for the constants.

    state_symbols = (y, uy)
    specified_symbols = (u,)

    # Specify the objective function and form the gradient.
    obj_func = sm.Integral(y**2 + u**2, t)
    sm.pprint(obj_func)
    obj, obj_grad = create_objective_function(obj_func,
                                          state_symbols,
                                          specified_symbols,
                                          tuple(),
                                          num_nodes,
                                          node_time_interval=interval_value)

    # Specify the symbolic instance constraints, as per the example
    instance_constraints = (
        y.func(t0) - 1,
        y.func(tf) - 1.5,
    )

    # Create the optimization problem and set any options.
    prob = Problem(obj, obj_grad, eom, state_symbols,
               num_nodes, interval_value,
               instance_constraints=instance_constraints,
    )

    prob.add_option('nlp_scaling_method', 'gradient-based')

    # Give some rough estimates for the trajectories.
    initial_guess = np.zeros(prob.num_free)

    # Find the optimal solution.
    solution, info = prob.solve(initial_guess)
    print(info['status_msg'])
    print(f'Objective value achieved: {info['obj_val']:.4f}, as per the book ' +
        f'it is {6.7241}, so the error is: '
        f'{(info['obj_val'] - 6.7241)/6.7241*100:.3f} % ')

    # Plot the optimal state and input trajectories.
    prob.plot_trajectories(solution)

    # Plot the constraint violations.
    prob.plot_constraint_violations(solution)

    # Plot the objective function as a function of optimizer iteration.
    prob.plot_objective_value()

# %%
# As per the example tf = 10000

tf = 10000
num_nodes = 501
solve_optimization(num_nodes, tf)

# %%
# With the value of tf = 10000 above, opty converged to a locally optimal point,
# but the objective value is far from the one given in the book.
# As per the plot of the solution y(t) it seems, that most of the time y(t) = 0,
# only at the very beginning and the very end it is different from 0.
# So, it may make sense to use a smaller tf.
# Also increasing num_nodes may help.

tf = 8.0
num_nodes = 10001
solve_optimization(num_nodes, tf)

# %%

# sphinx_gallery_thumbnail_number = 4

plt.show()

[moorepants]

You need to open a PR with the new example if you want help with it. Otherwise we can't run it in context or see what errors the CI produces.

[Peter230655]

Will do so tomorrow morning.

[Peter230655]

I opened PR 256 this morning.

[Peter230655]

Problem was the wrong file name. For examples-gallery it must start with plot_...

[moorepants]

Reopened, this is still an issue we can solve.

[Peter230655]

Reopened, this is still an issue we can solve.

I thought you set up examples-gallery like this on purpose.

[moorepants]

This issue is about opty needing support for one scalar differential equation. When we add that support, we can close this issue.

[Peter230655]

This issue is about opty needing support for one scalar differential equation. When we add that support, we can close this issue.

This is what I get for mixing two things in the same issue! Sorry!