Bug tf variable

[ocots]

It seems that this is not correct:

julia> @def ocp begin
        tf ∈ R, variable
        t ∈ [0, tf], time
        x ∈ R, state
        u ∈ R, control
        ẋ(t) == tf * u(t)
        x(0) == 0
        x(tf) == 1
        tf + 0.5∫(u(t)^2) → min
        end

The (autonomous) optimal control problem is given by:

    tf ∈ R, variable
    t ∈ [0, tf], time
    x ∈ R, state
    u ∈ R, control
    ẋ(t) == tf * u(t)
    x(0) == 0
    x(tf) == 1
    tf + 0.5 * ∫(u(t) ^ 2) → min

The (autonomous) optimal control problem is of the form:

    minimize  J(x, u, tf) = g(x(0), x(tf), tf) + ∫ f⁰(x(t), u(t), tf) dt, over [0, tf]

    subject to

        ẋ(t) = f(x(t), u(t), tf), t in [0, tf] a.e.,

        ϕl ≤ ϕ(x(0), x(tf), tf) ≤ ϕu, 

    where x(t) ∈ R, u(t) ∈ R and tf ∈ R.

Declarations (* required):

╭────────┬────────┬──────────┬──────────┬───────────┬────────────┬─────────────╮
│ times* │ state* │ control* │ variable │ dynamics* │ objective* │ constraints │
├────────┼────────┼──────────┼──────────┼───────────┼────────────┼─────────────┤
│   V    │   V    │    V     │    V     │     V     │     V      │      V      │
╰────────┴────────┴──────────┴──────────┴───────────┴────────────┴─────────────╯


julia> sol = solve(ocp)
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:     1003
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:      203

Total number of variables............................:      304
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:      203
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  2.0000000e-01 9.00e-01 4.91e-03   0.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  2.3319855e-01 8.99e-01 3.84e-01 -11.0 3.36e+01    -  1.00e+00 9.77e-04h 11
   2  3.1063422e-01 8.96e-01 1.15e+00 -11.0 1.94e+01    -  1.00e+00 3.91e-03h  9
   3  5.7872626e-01 8.82e-01 1.03e+00 -11.0 1.74e+01  -2.0 1.00e+00 1.56e-02h  7
   4  8.6337333e-01 8.68e-01 7.24e-01 -11.0 1.83e+01  -0.7 1.00e+00 1.56e-02h  7
   5  1.2396253e+00 8.14e-01 1.02e+00 -11.0 6.07e+00    -  1.00e+00 6.25e-02h  5
   6  8.3728076e-01 7.10e-03 4.56e+00 -11.0 9.76e-01  -1.1 1.00e+00 1.00e+00h  1
   7  1.6241923e+00 1.05e-02 6.93e+00 -11.0 1.43e+00    -  1.00e+00 1.00e+00h  1
   8  1.5839909e+00 2.95e-04 6.23e-01 -11.0 1.00e-01  -0.7 1.00e+00 1.00e+00h  1
   9  1.2672905e+00 1.12e-02 4.91e-01 -11.0 1.95e+00    -  1.00e+00 1.00e+00f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.3583706e+00 2.54e-03 3.59e-03 -11.0 6.11e-01    -  1.00e+00 1.00e+00h  1
  11  1.4778821e+00 1.13e-04 9.96e-03 -11.0 9.02e-02    -  1.00e+00 1.00e+00h  1
  12  1.4755675e+00 1.90e-07 1.54e-06 -11.0 4.38e-03    -  1.00e+00 1.00e+00h  1
  13  1.4755759e+00 5.11e-12 2.99e-10 -11.0 2.27e-05    -  1.00e+00 1.00e+00h  1

Number of Iterations....: 13

                                   (scaled)                 (unscaled)
Objective...............:   1.4755758928406426e+00    1.4755758928406426e+00
Dual infeasibility......:   2.9935254275414991e-10    2.9935254275414991e-10
Constraint violation....:   5.1116888499791457e-12    5.1116888499791457e-12
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   0.0000000000000000e+00    0.0000000000000000e+00
Overall NLP error.......:   2.9935254275414991e-10    2.9935254275414991e-10


Number of objective function evaluations             = 47
Number of objective gradient evaluations             = 14
Number of equality constraint evaluations            = 53
Number of inequality constraint evaluations          = 0
Number of equality constraint Jacobian evaluations   = 14
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations             = 13
Total seconds in IPOPT                               = 0.190

EXIT: Optimal Solution Found.
CTBase.OptimalControlSolution

julia> p0 = sol.costate(0)
0.7377879463735035

julia> tf = sol.variable
1.1066819197577107

[ocots]

With an indirect method I get:

julia> H(x, p, u, tf) = tf*p*u-u^2/2
H (generic function with 2 methods)

julia> u(x, p, tf) = tf*p
u (generic function with 1 method)

julia> H(x, p, tf) = H(x, p, u(x, p, tf), tf)
H (generic function with 2 methods)

julia> F = Flow(ocp, u);

julia> function S(s, p0, tf)
       xf, pf = F(0, 0, p0, tf, tf)
       s .= [xf-1; H(xf, pf, tf)-1]
       end
S (generic function with 2 methods)

julia> nle = (s, ξ, λ) -> S(s, ξ[1], ξ[2])
#19 (generic function with 1 method)

julia> prob = NonlinearProblem(nle, ξ)
NonlinearProblem with uType Vector{Int64}. In-place: true
u0: 2-element Vector{Int64}:
 1
 1

julia> indirect_sol = NonlinearSolve.solve(prob)
retcode: Success
u: 2-element Vector{Float64}:
 1.6817928305074286
 0.8408964152537142

julia> p0_sol = indirect_sol.u[1]
1.6817928305074286

julia> tf_sol = indirect_sol.u[2]
0.8408964152537142

[ocots]

It was a problem to make some tests to: control-toolbox/CTFlows.jl#15.

[jbcaillau]

well mr @ocots, I would first put the problem in canonical form (as such, it is not), having $t_f$ as an additional state variable, and writing the PMP (and the shooting function) properly 🙂

[PierreMartinon]

Hi Olivier, JB

This is actually a nice test problem with a Bolza objective. There are none currently in the test suite so bugs may have slipped through. Also, until very recently the Mayer objective was in autonomous form only so there may be something here.

I'll check this !
Pierre

[ocots]

well mr @ocots, I would first put the problem in canonical form (as such, it is not), having $t_f$ as an additional state variable, and writing the PMP (and the shooting function) properly 🙂

With PMP I got the same as with indirect shooting.

[PierreMartinon]

With PMP I got the same as with indirect shooting.

@ocots can you plot the indirect solution ?
The direct one:

Also I noticed we have the free tf explicitely in the dynamics, this is new as well.

[ocots]

Indirect solution

[PierreMartinon]

Thanks !
The objective for the direct solution is 1.48, consistent with the values of u=0.82 and tf=1.107.

For the indirect solution the objective seems to be 1.67 if we take tf=0.84 and u=1.41. Could we have a local solution ? Can you try a solve with the direct data as initial guess ?

Ps. @jbcaillau the non-autonomous Mayer is not yet in CTBase main right ?

[ocots]

If I am not wrong, the solution is $t_f = (1/2)^{1/4}$ and $p_0=2t_f$ in the case where $p^0=-1$. No local solution I think.

[jbcaillau]

using OptimalControl
using NLPModelsIpopt 
using OrdinaryDiffEq

@def ocp begin
    tf ∈ R, variable
    t ∈ [0, tf], time
    x ∈ R, state
    u ∈ R, control
    ẋ(t) == tf * u(t)
    x(0) == 0
    x(tf) == 1
    tf + 0.5∫(u(t)^2) → min
end

sol = solve(ocp)
x = sol.state
u = sol.control
tf = sol.variable
p0 = sol.costate(0)

@def ocp2 begin
    s ∈ [0, 1], time
    y = (x, tf) ∈ R², state
    u ∈ R, control
    dx = tf(s)^2 * u(s)
    dtf = 0 * u(s) # 0
    ẏ(s) == [dx, dtf]
    x(0) == 0
    x(1) == 1
    tf(1) + 0.5∫(tf(s) * u(s)^2) → min
end

init2 = (state=s -> [x(tf * s), tf], control=s -> u(tf * s))
sol2 = solve(ocp2; init=init2)
tf2 = sol2.state(0)[2]
p02 = sol2.costate(0)[1]

ur(y, q) = y[2] * q[1]
f = Flow(ocp2, ur)

function S(tf, p0)
    yf, qf = f(0, [0, tf], [p0, 0], 1)
    s = [yf[1] - 1, qf[2] + 1]
    return s
end

S(tf2, p02)

[jbcaillau]

writing the PMP (and the shooting function) properly 🙂

[ocots]

@jbcaillau Alors, tu trouves quoi ?

Après, moi je ne veux pas transformer le problème car j'avais besoin d'un exemple pour ça : control-toolbox/CTFlows.jl#38

[jbcaillau]

Essaie tu verras 🫢

[jbcaillau]

@PierreMartinon #190

[PierreMartinon]

If I am not wrong, the solution is tf=(1/2)1/4 and p0=2tf in the case where p0=−1. No local solution I think.

@jbcaillau Alors, tu trouves quoi ?

Après, moi je ne veux pas transformer le problème car j'avais besoin d'un exemple pour ça : control-toolbox/CTFlows.jl#38

I tried to check the dynamics and final state for both our solutions, and both seem correct oO
By the way this is a nice test problem indeed !

[ocots]

Essayez de mettre 200 pas.

[PierreMartinon]

Ok vous m'avez perdu :D
Concretement, on en est ou ?

(avec 200 steps je retrouve la meme solution qu'a 100, par contre il faut ajouter une borne inf sur tf, comme a chaque fois qu'on est a tf libre en fait, sinon on risque toujours d'aller dans les tf negatifs)

[ocots]

En effet, en utilisant proprement le PMP on retrouve bien le bon résultat :-)

La fonction de tir :

function S(s, p0, tf)
    xf, pf = F(0, 0, p0, tf, tf)
    s .= [xf-1; H(xf, pf, tf)+tf^2*pf^2-1]
end

et on obtient :

p0 =  0.7377879464668851
tf =  1.1066819197003364

A la main, la solution est : $t_f = (3/2)^{1/4}$ et $p_0 = 2t_f/3$.

[ocots]

@PierreMartinon Ici tu as 3 exemples avec t0, tf et les deux libres. Avec les temps qui interviennent dans la dynamique, le coût, etc : https://github.com/control-toolbox/CTFlows.jl/pull/38/files#diff-6245c325415158ae48d711d73ce595e7d2dfb96f7bf04341fe10ef35952c8bf1.

[jbcaillau]

@PierreMartinon #190

@PierreMartinon please take a look

[PierreMartinon]

@PierreMartinon Ici tu as 3 exemples avec t0, tf et les deux libres. Avec les temps qui interviennent dans la dynamique, le coût, etc : control-toolbox/CTFlows.jl#38 (files).

Pssst, celui avec t0,tf libres a l'air d'admettre une autre (meilleure) solution avec t0=0.14 :D

[ocots]

@PierreMartinon Ici tu as 3 exemples avec t0, tf et les deux libres. Avec les temps qui interviennent dans la dynamique, le coût, etc : control-toolbox/CTFlows.jl#38 (files).

Pssst, celui avec t0,tf libres a l'air d'admettre une autre (meilleure) solution avec t0=0.14 :D

J'avoue que je ne l'ai pas résolu même si j'ai écrit solution.