Error message

[Peter230655]

I am trying to model a differential game. (War of Attrition and Pursuit)
I have to set up the objective function and its gradient manually.
Setting up Problem seems fine, but when I run solve, I get this error message:

Would this indicate a mistake in my gradient of the objective function?
My objective function depends on a known_trajectory_map, and of course so does its gradient.
Is this allowed?
Thanks!

[moorepants]

My first guess is that your EoMs or their Jacobian returns invalid values. You can manually check whether those evaluate before you call solve().

[moorepants]

See the documentation for how to evaluate the Jacobian of the constriants: https://opty.readthedocs.io/stable/api.html#opty.direct_collocation.Problem

[Peter230655]

See the documentation for how to evaluate the Jacobian of the constriants: https://opty.readthedocs.io/stable/api.html#opty.direct_collocation.Problem

I am so hopelessly dumb sometimes! :-(
I thought I need this

and I thought it was a method of ConstraintCollocator, but it is not.
How do I use it?
Thanks!

[moorepants]

You don't need to generate a Jacobian you just evaluate it. See the jacobian() method.

[Peter230655]

You don't need to generate a Jacobian you just evaluate it. See the jacobian() method.

Now I got you - age!!
I did this:

eom = sm.Matrix([ux1 - (m1 - c1 * phsi(t) * x2), ux2 - (m2 - c2 * phi(t) * x1)])
instance_constraints = (
    x1.subs({t: t0}) - 10,
    x2.subs({t: t0}) - 10,
)
constraints = eom.col_join(sm.Matrix(instance_constraints))
jakob = constraints.jacobian([x1, x2 ,ux1, ux2] + [phi(t), phsi(t)])

and this result came:

It seems to give $\dfrac{d}{d \phi (t)} (\phi(t0) = 0$ which is correct, I believe, since $\phi(t0) = constant$

[moorepants]

Now I got you - age!!

Sorry, you haven't got it. Just call Problem.jacobian(values) to see what the Jacobian evaluation returns.

[Peter230655]

I am 69 1/2 - by the time I'll be 75 I may be o.k....
I did this:

jakob = problem_phi.jacobian(initial_guess_phi)
print('max entry in jakob', np.max(jakob))
print('min entry in jakob', np.min(jakob))
print('len(jakob)', len(jakob))

and got these resulst;

max entry in jakob 3.266960943401931
min entry in jakob -2.8098129349783894
len(jakob) 884

does this look o.k.?

[moorepants]

I can't tell because you only show two values. If there are ever any nans or infs in that it will fail.

[Peter230655]

I can't tell because you only show two values. If there are ever any nans or infs in that it will fail.

I printed out the complete jakob and saw no nans of infs.
I do not see how the length of 884 comes up:

I set num_nodes = 50
there are two oems and two instance_constraints.
There are six variables (x1, x2, ux1, ux2, phi, phsi). phsi is given as a known_trajectory_map

No matter how I divide it, I do not get an integer.

[moorepants]

It is difficult to help without sample code that I can run.

[Peter230655]

Let me clean up the code and send it tomorrow morning.
The basic idea is this:
two players are playing a differential game.
Player 1 has strategy phi, player 2 has strategy phsi.
There are oems depending on phi and phsi.
There is a payoff function P which player 1 wants to minimize and player 2 wants to maximise.
P depends on phi, phsi and on the solution of the oems.

So, I set up two Problems:
problem_phi uses the strategy phsi as a known_trajectory and tries minimize P
problem phsi uses the strategy phi as known_trajectory and tries to minimize -P

My hope is, that by iterating this it will converge to optimal strategies phi and phsi.
(My book says, that they exist. They are discontinuous, but I am not bothered at the moment.)

In which form would you like the code?
Just in .....in here?
Thanks for taking so much time!

[Peter230655]

Below is the program.
Thanks again for taking the time!!

# %%
"""
War of Attrition and Attack
===========================

Two players, called :math:`\phi` and :math:`\psi`, have weapons of quantity
$:math`x_1` and :math:`x_2`, respectively. They may use them to either destroy
the other player's weapons or to attack the other player. They employ strategfies
:math:`\phi(t)` and :math:`\psi(t)` with :maht:`0 \leq \phi(t), \psi(t) \leq 1`.
The amount of weapons left is governed by:

 :math:`\dot{x}_1 = m_1 - c_1 \psi(t) x_2`
 :math:`\dot{x}_2 = m_2 - c_2 \phi(t) x_1`

The pay off is:

 :math:` P(x, \psi, \phi) = \int_0^T (1 - \psi(t)) x_2 - (1 - \phi(t)) x_1 dt`

Player :math:`\phi` wants to maximize the payoff and player :math:`\psi` wants to
minimize it.

The idea comes from the book Differential Games by Avner Friedman, page 26.

The idea is this: player :math:`\phi` wants to minimize P(..), so this will be
his objective function and the oppent's strategie :math:`psi` enters as a
known_trajectory_map.  Player :math:`\psi` wants to maximize P(..), so his
objective function will be -P(..) and the opponent's strategie :math:`\phi`
enters as a known_trajectory_map.

By iterating a few times maybe the strategies will converge to a solution.
/AS per the book, there is a solution, but it is discontinuous

"""
# %%
# Here I simply integrate the eoms, just to see what they may look like
#
# %%
import sympy as sm
import numpy as np
import sympy.physics.mechanics as me
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
print('Sympy version: ', sm.__version__)
from opty import Problem

m1, m2, c1, c2, psi, phi = sm.symbols('m1 m2 c1 c2 psi phi')

par_map = {
    m1: 1,
    m2: 1,
    c1: 1.,
    c2: 1.,
    psi: 0.25,
    phi: 0.255,
}

def gradient(t, x, args):
    A = args[0] - args[2] * args[4] * x[1]
    B = args[1] - args[3] * args[5] * x[0]
    return np.array([A, B])

t0, tf = 0, 10
num_nodes = 100
times = np.linspace(t0, tf, num_nodes)
t_span = (t0, tf)

x0 = np.array([10, 10])
pL_vals = [par_map[i] for i in par_map.keys()]

resultat1 = solve_ivp(gradient, t_span, x0, t_eval = times, args=(pL_vals,)) #, method='Radau')
resultat = resultat1.y.T
print('Shape of result: ', resultat.shape)
event_dict = {-1: 'Integration failed', 0: 'Integration finished successfully', 1: 'some termination event'}
print(event_dict[resultat1.status], ' the message is: ', resultat1.message)

fig, ax = plt.subplots()
plt.plot(times, resultat[:, 0], label='x1')
plt.plot(times, resultat[:, 1], label='x2')
plt.xlabel('Time')
ax.legend();

plt.show()

# %%
# Set up the Optimization Problem.
# --------------------------------
#
x1, x2, ux1, ux2 = me.dynamicsymbols('x1 x2 ux1 ux2')#
psi, phi = sm.symbols('psi phi', cls=sm.Function)
m1, m2, c1, c2 = sm.symbols('m1 m2 c1 c2')

t = me.dynamicsymbols._t
t0, tf = 0, 10
num_nodes = 50

eom = sm.Matrix([ux1 - (m1 - c1 * psi(t) * x2), ux2 - (m2 - c2 * phi(t) * x1)])
instance_constraints = (
    x1.subs({t: t0}) - 10,
    x2.subs({t: t0}) - 10,
)
constraints = eom.col_join(sm.Matrix(instance_constraints))
jakob = constraints.jacobian([x1, x2 ,ux1, ux2] + [phi(t), psi(t)])

par_map = {
    m1: 2,
    m2: 2,
    c1: 1.,
    c2: 1,
}


state_symbols = [x1, x2, ux1, ux2]
constant_symbols = [m1, m2, c1, c2]
interval_value = (tf - t0) / (num_nodes - 1)

"""
objective_phi = sm.Integral(((1 - psi(t))*x2 - (1 - phi(t))*x1), t)
objective_psi = -objective_phi
"""

phi_dict = {phi(t): np.ones(num_nodes)*0.25}
psi_dict = {psi(t): np.ones(num_nodes)*0.25}
initial_guess_phi = np.random.randn(5*num_nodes) * 1
initial_guess_psi = np.random.randn(5*num_nodes) * 1

def obj_phi(free):
    h1 = -interval_value * np.sum([(1 - free[4*num_nodes+i])*free[i] for i in range(num_nodes)])
    h2 = interval_value * np.sum([(1 - psi_dict[psi(t)][i])*free[num_nodes+i] for i in range(num_nodes)])
    return h1 + h2

def obj_psi(free):
    h1 = interval_value * np.sum([(1 - phi_dict[phi(t)][i])*free[i] for i in range(num_nodes)])
    h2 = -interval_value * np.sum([(1 - free[4*num_nodes+i])*free[num_nodes+i] for i in range(num_nodes)])
    return h1 + h2

def obj_grad_phi(free):
    grad = np.zeros_like(free)
    grad[:num_nodes] = -interval_value * (1 - free[4*num_nodes:])
    grad[num_nodes:2*num_nodes] = interval_value * (1 - psi_dict[psi(t)])
    return grad

def obj_grad_psi(free):
    grad = np.zeros_like(free)
    grad[: num_nodes] = interval_value * (1 - phi_dict[phi(t)])
    grad[num_nodes:2*num_nodes] = -interval_value * (1 - free[4*num_nodes:])
    return grad

instance_constraints = (
    x1.subs({t: t0}) - 10,
    x2.subs({t: t0}) - 10,
)
constraints = eom.col_join(sm.Matrix(instance_constraints))
jakob = constraints.jacobian(state_symbols + [phi(t), psi(t)])

bound_phi = {phi(t): (0, 1)}
bound_psi = {psi(t): (0, 1)}

problem_psi = Problem(
    obj_psi,
    obj_grad_psi,
    eom,
    state_symbols,
    num_nodes,
    interval_value,
    known_parameter_map=par_map,
    known_trajectory_map=phi_dict,
    instance_constraints=instance_constraints,
    bounds=bound_psi,
    time_symbol=me.dynamicsymbols._t,
)
print('psi finished')

problem_phi = Problem(
    obj_phi,
    obj_grad_phi,
    eom,
    state_symbols,
    num_nodes,
    interval_value,
    known_parameter_map=par_map,
    known_trajectory_map=psi_dict,
    instance_constraints=instance_constraints,
    bounds=bound_phi,
    time_symbol=me.dynamicsymbols._t,
)
print('phi finished')

jakob = problem_phi.jacobian(initial_guess_phi)
print('max entry in jakob', np.max(jakob))
print('min entry in jakob', np.min(jakob))
print('len(jakob)', len(jakob))
print('jakob:', '\n')
print(jakob)

# %%
# Solve the Optimization Problem.
# -------------------------------

for i in range(1):

    solution_phi, info = problem_phi.solve(initial_guess_phi)
    print(info['status_msg'])
    problem_phi.plot_objective_value()
    initial_guess_phi = solution_phi
    phi_dict = {phi(t): solution_phi[4*num_nodes :]}
    print(i, 'phi')

    solution_psi, info = problem_psi.solve(initial_guess_psi)
    print(info['status_msg'])
    problem_psi.plot_objective_value
    psi_dict = {psi(t): solution_psi[4*num_nodes :]}
    initial_guess_psi = solution_psi
    print(i, 'psi')

[Peter230655]

Another idea is to try to solve all of the examples in Betts' book. I think opty can solve some of them but we'd find many weaknesses due to lack of features for some of them. Pycollo also has a bunch of example problems: https://github.com/brocksam/pycollo/tree/main/examples most of which opty cannot solve. You could play around with pycollo too.

I had shyed away from Betts' book as it cost close to € 400.-
Just now I found the latest (2020) edition at € 130.- so I bought it.
Of course I read about pycollo, but not really looked at it. I will do this.

[moorepants]

I suspect pycollo can solve all problem's in Bett's book. It was more-or-less designed to implement Bett's methods. opty only implements the most basic aspects of Bett's methods. pycollo has no documentation though and it chokes on very large EoMs (like the bicycle model).

[Peter230655]

I suspect pycollo can solve all problem's in Bett's book. It was more-or-less designed to implement Bett's methods. opty only implements the most basic aspects of Bett's methods. pycollo has no documentation though and it chokes on very large EoMs (like the bicycle model).

If you think, how helpless I am with good documentation, pycollo might not be for me.....
Betts' book will arrive in two weeks, and until then I will be busy with the two suggestios you gave me. :-)
NB: it is now running fine, but does not seem to converge.

[moorepants]

If you think, how helpless I am with good documentation, pycollo might not be for me.....

I expected that :)

[Peter230655]

If you think, how helpless I am with good documentation, pycollo might not be for me.....

I expected that :)

:-) :-) :-)
I will certainly look at it!

[Peter230655]

Enough to keep me going for a while!
THANKS!!

[moorepants]

You can open an issue about this to request a better error message if the user does not pass in differential algebraic equations. We could have some kind of input check.

[Peter230655]

You can open an issue about this to request a better error message if the user does not pass in differential algebraic equations. We could have some kind of input check.

I will do this.
If I understand correctly, the equations_of_motion must contain differential equations, too, no only AEs?
This test would be in the class Problem?

[moorepants]

I think they have to be DAEs or ODEs and if you only have AEs, then it probably fails. If you want to solve only AEs, then I don't think the direct collocation method is needed.

[Peter230655]

I may have found a simple way of checking:

eom = sm.Matrix([a.diff(t) - b.diff(t,2) + c.diff(t), a.diff(t, 2)])
for a in sm.preorder_traversal(eom):
    a = str(a)
    if 'Derivative' in a:
        print(a, 'there is a derivative')

almost too simple to be true...
I will play around wqith it a bit more.

[Peter230655]

expr.has(sm.Derivative) should be even faster.

MUCH faster! I cannot check the time anymore.

[tjstienstra]

It does use caching, so it may cheat with your time results:

expr.has(sm.Derivative)   # Takes a short moment
expr.has(sm.Derivative)  # Is instant, as the value comes from the cache
expr.has.cache_clear()
expr.has(sm.Derivative)   # Now it actually has to do something ;)

[Peter230655]

It does use caching, so it may cheat with your time results:

expr.has(sm.Derivative)   # Takes a short moment
expr.has(sm.Derivative)  # Is instant, as the value comes from the cache
expr.has.cache_clear()
expr.has(sm.Derivative)   # Now it actually has to do something ;)

expr.has.cache_clear()
raises an error in my Jupyter Notebook:
AttributeError: 'function' object has no attribute 'cache_clear'

Background: I tried an opty simulation and through stupidity I only gave algebraic equations as eoms. This raised an error which was strange. When Jason found it, he suggested one should test for this, and give a suitable error message.

@moorepants: Should I make a PR to add Timo's line? It seems lighteningly fast. If YES, it would be neaqr the top of the class Problem?

[tjstienstra]

expr.has.cache_clear()
raises an error in my Jupyter Notebook:
AttributeError: 'function' object has no attribute 'cache_clear'

Hmmm, thinks give a bit different results then I expected (something doesn't seem to work entirely like I expected):

from timeit import timeit

import sympy as sm

a, b, c, x = sm.symbols("a b c x")
f = sm.Function("f")(x)
expr = ((a * b + f.diff())**20).expand() * (c + f)
print(expr.has(sm.Derivative))
n = 100000
print(timeit(lambda: expr.has(sm.Derivative), number=n))  # 0.02411990000109654
print(expr.has.cache_info())  # CacheInfo(hits=100005, misses=46, maxsize=1000, currsize=46)
print(timeit(lambda: expr.has(sm.Derivative), lambda: expr.has.cache_clear(), number=n))  # 0.023935599998367252
print(expr.has.cache_info())  # CacheInfo(hits=99999, misses=1, maxsize=1000, currsize=1)

Background: I tried an opty simulation and through stupidity I only gave algebraic equations as eoms. This raised an error which was strange. When Jason found it, he suggested one should test for this, and give a suitable error message.

Thanks for the explanation.

[Peter230655]

In any case lighteningly fast, compared to my clumsy idea!! :-)

[Peter230655]

A check for algebraic equations of motion was added to opty.