Function to generate object and objective gradient evaluator functions

[tjstienstra]

This PR proposes a utility function to allow one to specify an analytic expression as the objective function. It is a relatively simple function with one major limitation, namely that one cannot use time-dependent and time-independent variables in the objective function.

Additionally, I noticed while opening this PR that #31 aims to implement a similar utility in an entirely different way. Will have to look into that one to see if it would be a better solution.

[tjstienstra]

Okay, based on the other PR. I would say that this is a fun function for quite a few cases. However, we optimally want to define the objective expression as a continuous function, possibly with an integral.

import sympy as sm
import sympy.physics.mechanics as me

x, v, f = me.dynamicsymbols("x v f")
t = me.dynamicsymbols._t
m = sm.symbols("m")
integration_method = "backward euler"
node_times = np.linspace(0, 2, 101)  # duration=2.0 ; N=101 ; interval_value=0.02

obj_expr = m**2 + sm.Integral(f**2, (t, 0.5, 1.5))
obj, obj_grad = create_objective_function(obj_expr, (s, v), (f,), (m,), integration_method, node_times)

# Should give something like
# free[303] is just the value of m
# free[228] is f.subs(t, 0.52) (we are using backward euler), therefore it is not f.subs(t, 0.5)
# free[278] is f.subs(t, 1.50)
def obj(free):
    return free[303] ** 2 + np.sum(free[228:278] ** 2)

def obj_grad(free):
    return np.hstack((np.zeros(101), np.zeros(101), np.zeros(25), 2 * free[228:278], np.zeros(26), 2 * free[303]))

Not yet fully sure if this is the ideal method to construct obj_grad, but it is the version which looks most like what this PR currently tried to do.

[tjstienstra]

A temporary simplification would be to just ignore the bounds of the integral for now, i.e. assert integral.limits == ((t,),)

[tjstienstra]

Here is a fun objective, that will give even when trying to keep things simple quite a few problems.

import numpy as np
import sympy as sm
import sympy.physics.mechanics as me

x, v, f = me.dynamicsymbols("x v f")
t = me.dynamicsymbols._t
m, k, c = sm.symbols("m k c")
integration_method = "backward euler"
node_times = np.linspace(0, 2, 11)

states = sm.ImmutableDenseMatrix([x, v])
inputs = sm.ImmutableMatrix([f])
params = sm.ImmutableMatrix([c, k, m])
n, q, r = states.shape[0], inputs.shape[0], params.shape[0]

obj = m**2 + sm.Integral(f**2, (t,)) + sm.Integral(c**2 * f, (t,)) + sm.Integral(k**2 , (t,))
obj_grad = sm.ImmutableMatrix([obj]).jacobian(states.col_join(inputs).col_join(params))

[tjstienstra]

@moorepants I think I figured out a method to combine it all. It is not the most beautiful function I have ever written, but it seems to do the job as the tests show.

[moorepants]

https://github.com/csu-hmc/opty/pull/116 fixes the test failuers (pytest 8 just came out with a deprecation). You can merge master.

[moorepants]

You happy with this? I think it is good to go in if you do.

[tjstienstra]

I think that it is good for now. Unless you're planning to add a lot of advanced algorithms soon, I really don't see the need to make it more advanced. Like the internals may change. The only thing is that the keyword arguments expect fixed time stepping. All other limitations seem to be covered using errors.

[moorepants]

Ok. let's merge. Nice addition!