Anciens codes

[ocots]

J'avais écrit des codes pour le problème de saturation et le problème de contraste. Je vous les joins juste après.

Remarque. Ce serait bien de pouvoir essayer avec grille non uniforme pour améliorer la précision du second bang dans le problème de saturation. Voici une figure pour le problème de saturation (cas du sang) avec N = 50 puis N = 500 points.

Pour le problème de contraste, j'obtiens un contraste à 0.474 pour un temps à $1.5 T_\min$ :

[ocots]

Problème de saturation :

using OptimalControl
using Plots

# Blood case
Tmin = 2π*6.7981
T₁ = 1.35 # s
T₂ = 0.05

@def ocp begin

    tf ∈ R, variable
    t ∈ [ 0, tf ], time
    q = [ y, z ] ∈ R², state
    u ∈ R, control

    tf ≥ 0
    q(0)  == [ 0, 1 ]
    q(tf) == [ 0, 0 ]

    -1 ≤ u(t) ≤ 1

    q̇(t) == F₀(q(t)) + u(t) * F₁(q(t))

    tf → min

end

#
ω = 2π*32.3 # Hz
γ = 1/(ω*T₁)
Γ = 1/(ω*T₂)

#
F₀(q) = begin
    y, z = q
    [ -Γ * y, γ * (1 - z) ]
end

F₁(q) = begin
    y, z = q
    [ -z, y ]
end

# plot
function spinplot(sol, N)

    tf = sol.variable
    y = t -> sol.state(t)[1]
    z = t -> sol.state(t)[2]
    u = sol.control

    tspan  = range(0, tf, N)
    yz_plt = plot(y.(tspan), z.(tspan), xlims=(-1, 1), ylims=(-1, 1), xlabel="y", ylabel="z", label=:none)
    u_plt  = plot(tspan, u.(tspan), ylims=(-1.1, 1.1), xlabel="t", ylabel="u", label=:none)
    return plot(yz_plt, u_plt, layout=(1, 2), size=(800, 400))

end

# solve
N = 50
sol = solve(ocp, 
            init=(state=[ -0.9, 0.1 ], control=[ 0.05 ], variable=Tmin),
            grid_size=N, 
            tol=1e-6)

p1 = spinplot(sol, N)

#
N = 500
sol = solve(ocp, 
            init=sol,
            grid_size=N,
            tol=1e-12)

p2 = spinplot(sol, N)

# all plots
plot(p1, p2, layout=(2, 1), size=(800, 800))

[ocots]

Problème de contraste :

using OptimalControl
using Plots

function Constrast(F₀, F₁, tf)
    @def ocp begin

        t ∈ [ 0, tf ], time
        q = [ y₁, z₁, y₂, z₂ ] ∈ R⁴, state
        u ∈ R, control

        q(0)   == [ 0, 1, 0, 1 ]
        y₁(tf) == 0
        z₁(tf) == 0

        -1 ≤ u(t) ≤ 1

        q̇(t) == F₀(q(t)) + u(t) * F₁(q(t))

        ( y₂(tf)^2 + z₂(tf)^2 ) → max

    end
    return ocp
end

# Blood case
Tmin = 2π*6.7981
T₁_d = 1.35 # s
T₂_d = 0.05
T₁_o = 1.35 # s
T₂_o = 0.2

# Fluide case
#Tmin = 2π*26.1704
#T₁_d = 2.0 # s
#T₂_d = 0.3
#T₁_o = 2.5 # s
#T₂_o = 2.5

#
tf = 1.5*Tmin

#
ω    = 2π*32.3 # Hz
γ_d  = 1/(ω*T₁_d)
Γ_d  = 1/(ω*T₂_d)
γ_o  = 1/(ω*T₁_o)
Γ_o  = 1/(ω*T₂_o)

#
F₀(q) = begin
    y₁, z₁, y₂, z₂ = q
    [ -Γ_d * y₁, γ_d * (1 - z₁), -Γ_o * y₂, γ_o * (1 - z₂) ]
end

F₁(q) = begin
    y₁, z₁, y₂, z₂ = q
    [ -z₁, y₁, -z₂, y₂ ]
end

# plot
function spinplot(sol, N)

    q = sol.state
    y₁ = t -> q(t)[1]
    z₁ = t -> q(t)[2]
    y₂ = t -> q(t)[3]
    z₂ = t -> q(t)[4]

    u = t -> sol.control(t)

    tspan = range(0, tf, N)

    yz_1_plt = plot(y₁.(tspan), z₁.(tspan), xlims=(-1, 1), ylims=(-1, 1), xlabel="y₁", ylabel="z₁", label=:none)
    yz_2_plt = plot(y₂.(tspan), z₂.(tspan), xlims=(-1, 1), ylims=(-1, 1), xlabel="y₂", ylabel="z₂", label=:none)
    yz_plt   = plot(yz_1_plt, yz_2_plt, layout=(1, 2))
    u_plt  = plot(tspan, u.(tspan), ylims=(-1.1, 1.1), xlabel="t", ylabel="u", label=:none)

    return plot(yz_plt, u_plt, layout=(2, 1), size=(600, 600))

end

# ocp
ocp = Constrast(F₀, F₁, tf)

# solve
N = 100
sol = solve(ocp, 
            init=(state=[ 0, 0.2, 0, 0.8 ], control=0.05),
            grid_size=N, 
            tol=1e-12)

spinplot(sol, N)

# solve
N = 1000
sol = solve(ocp, 
            init=sol,
            grid_size=N, 
            tol=1e-12)

spinplot(sol, N)

[jbcaillau]

merci @ocots top (structure 3BS pour ce tf, à ce que je vois)

[ocots]

3 BS pour le contraste.

Pour le problème de saturation, ce serait bien d'avoir une grille plus fine au niveau du second arc bang qui est très court. Je n'arrive pas à bien le "capturer".