Bug with init from sol

[ocots]

There is bug in this example: https://github.com/control-toolbox/HMP/blob/main/main.ipynb

If you try to initialize with a previous solution.

See this part:

ρs = [0.1, 5, 10, 30, 100]
sols = []  
init = (state=t->[t, 1+t], control=[-1, 1], variable=0.5)
for ρ in ρs
    ocp = myocp(ρ)
    sol = solve(ocp, init=init)
    # init = sol # this is not working
    push!(sols, sol)
end

[PierreMartinon]

@ocots Shouldn't it be as below ? I'll try to add your example in CTDirect in test_continuation.jl

init = OptimalControlInit(sol)

[PierreMartinon]

@ocots Yup, it works. By the way it would be nice indeed not to have to call OptimalControlInit explicitely, maybe I can make it internal in solve ?

Te paramtric ocp

# definition of the parametric OCP
relu(x) = max(0, x)
μ = 10
p_relu(x) = log(abs(1 + exp(μ*x)))/μ
f(x) = 1-x
m(x) = (p_relu∘f)(x)
T = 2

function myocp(ρ)
    @def ocp begin
        τ ∈ R, variable
        s ∈ [ 0, 1 ], time
        x ∈ R², state
        u ∈ R², control
        x₁(0) == 0
        x₂(0) == 1
        x₁(1) == 1
        ẋ(s) == [τ*(u₁(s)+2), (T-τ)*u₂(s)]
        -1 ≤ u₁(s) ≤ 1
        -1 ≤ u₂(s) ≤ 1
        0 ≤ τ ≤ T
        -(x₂(1)-2)^3 - ∫( ρ * ( τ*m(x₁(s))^2 + (T-τ)*m(x₂(s))^2 ) ) → min
    end
    return ocp
end

The continuation

init3 = OptimalControlInit()
ρs = [0.1, 5, 10, 30, 100]
for ρ in ρs
    ocp3 = myocp(ρ)
    sol3 = solve(ocp3, print_level=0, init=init3)
    global init3 = OptimalControlInit(sol3)
    println("Rho ", ρ, " objective ", sol3.objective)
end

Result

Rho 0.1 objective -0.033579391539475464
Rho 5.0 objective -1.6863532482251995
Rho 10.0 objective -6.134496277986791
Rho 30.0 objective -34.55448970764526
Rho 100.0 objective -148.0223671431889

[ocots]

Yes it should be internal but it is nice to have the possibility to fo it externally.

[PierreMartinon]

Todo: avoid code duplication for the initial guess between OptimalControl and CTDirect
See what can be defined at the upper level (CTBase or OptimalControl) and then reused in CTDirect (and CTFlows).

@ocots I like your idea of having an upper level function that builds a generic functional initialization from the various possible inputs (vectors, functions, solution), used then by each method to generate its own specific initialization.

OCPInit(state=..., control=..., variable=...)

It would fuse the code currently in https://github.com/control-toolbox/OptimalControl.jl/blob/a7d6aa74d0cc702348861e8a3fb1e17f3e7bf2b4/src/solve.jl and https://github.com/control-toolbox/CTDirect.jl/blob/ed8c6e07126a08190ec53e1a05d4134f4c8fa466/src/init.jl , and I suggest we extend it to the costate as well (all arguments are supposed to be optional with a default constant initialization of 0.1 for instance). I would rather put OCPInit() in CTBase than OptimalControl since adding OptimalControl also adds CTFlows that is not really related to the direct method. A new source file init.jl ?

The actual initialization for the direct method is already in a second function CTDirect.initial_guess and would stay more or less the same. I suspect it would require only minimal changes in the indirect part as well ?

[ocots]

Just note that the init I use in OptimalControl.jl comes from CTDirect.jl

[PierreMartinon]

Yes I saw, this is why I said 'fuse' :D

[PierreMartinon]

Ok, I'll checkout CTBase on my side and try a first version with a PR (unless you want to do it ?)

[ocots]

You can do it :-)

Don't you think that this piece of code inside OptimalControl.jl/src/solve.jl should be placed directly inside the init function:

    if isnothing(init)
        init =  OptimalControlInit()
    elseif init isa CTBase.OptimalControlSolution
        init = OptimalControlInit(init)
    else
        x_init = :state    ∈ keys(init) ? init[:state]    : nothing
        u_init = :control  ∈ keys(init) ? init[:control]  : nothing
        v_init = :variable ∈ keys(init) ? init[:variable] : nothing
        init = OptimalControlInit(x_init=x_init, u_init=u_init, v_init=v_init)
    end

Of course, in this case, the function has to be updated.

[PierreMartinon]

Yes ! Ideally OCPInit will accept as many inputs as possible, and handle them to generate the functions for the variables. Then the module-specific functions for generating the actual initial guess and solving the problem can work in a simpler way. This will avoid duplications like when I started to handle in CTDirect the case of passing directly a solution as initialization, while you had already done it in CTBase.

[PierreMartinon]

First version is fine on the CTDirect side, I just replaced the former CTDirect.OptimalControlInit with a new CTBase.OCPInit.

mutable struct OCPInit

    state_init::Function
    control_init::Function
    variable_init::Union{Nothing, ctVector}
    costate_init::Function
    multipliers_init::Union{Nothing, ctVector}
    info::Symbol

    # constructor from constant vector or function handles
    function OCPInit(; state::Union{Nothing, ctVector, Function}=nothing, control::Union{Nothing, ctVector, Function}=nothing, variable::Union{Nothing, ctVector}=nothing, costate::Union{Nothing, ctVector, Function}=nothing, multipliers::Union{Nothing, ctVector}=nothing)

        init = new()
        init.info = :constant_or_function
        
        # use provided function or interpolate provided vector
        init.state_init = (state isa Function) ? t -> state(t) : t -> state
        init.control_init = (control isa Function) ? t -> control(t) : t -> control
        init.variable_init = variable
        
        init.costate_init = (costate isa Function) ? t -> costate(t) : t -> costate
        init.multipliers_init = multipliers

        return init
    end

    # constructor from existing solution
    function OCPInit(sol::OptimalControlSolution)

        init = new()
        init.info = :solution

        # Notes
        # - the possible internal state for Lagrange cost is not taken into account here
        # - set scalar format for dimension 1 case
        init.state_init    = t -> sol.state(t)
        init.control_init  = t -> sol.control(t)
        init.variable_init = sol.variable
        #+++ add costate and scalar multipliers

        return init
    end
end

@ocots about the nice code you mentioned, we could simply check in solve if the init passed is a OCPInit struct, and if not call something like `init=OCPInit(init)*. Then we can pass to solve either an OCPInit struct, or whatever the constructor OCPInit() accepts. I'll try that.

[PierreMartinon]

Ok, seems to work. When calling solve we can use both the explicit call to OCPInit or just pass the data for the initial guess.

For instance, a mixed functional x / constant u / default v initial guess

sol = solve(ocp, init=OCPInit(state=x_func, control=u_const))
sol = solve(ocp, init=(state=x_func, control=u_const))

and similarly for a warm start

sol = solve(ocp, init=OCPInit(sol0))
sol = solve(ocp, init=sol0)

Same for the other functions that take an initialization (directTrancription, setDOCPInit and solve applied to a DOCP instead of a OCP)

Changes are in the 'init_CT...' branches for CTBase and CTDirect. I think the only changes left to do are in OptimalControl. @ocots, do you want to do it ? In the end we'll need to do the merge and release for all 3 packages.

[ocots]

A first release of CTDirect must be done updating the Y of X.Y.Z version. Then, it will be simple to uodate OptimalControl.jl. I will do the update for OptimalControl.jl.

[PierreMartinon]

Release done for CTBase 0.8.0 and CTDirect 0.5.0

[ocots]

OK I do asap the release of OptimalControl.jl

[ocots]

Release done but I forgot to export new functions from CTDirect that can be useful.

[PierreMartinon]

We'll do it next time :-)