|
| 1 | +using LinearAlgebra |
| 2 | +using Statistics |
| 3 | + |
| 4 | +θ = [0.1,GaussNum(0.2,0.1),0.3,0.2] |
| 5 | +u0 = [GaussNum(1.0,0.1),GaussNum(1.0,0.1)] |
| 6 | +tspan = (0.,100.) |
| 7 | +dt = 0.1 |
| 8 | +prob = ODEProblem(f,tspan,u0,θ) |
| 9 | + |
| 10 | + |
| 11 | +isuncertain(x::GaussNum) = x.σ!=0 |
| 12 | +isuncertain(x::Number) = false |
| 13 | +sig(x::GaussNum) = x.σ |
| 14 | +mu(x::GaussNum) = x.μ |
| 15 | + |
| 16 | +#struct UncertainODEProblem{OP<:ODEProblem,S0<:AbstractMatrix,S<:AbstractMatrix,X<:AbstractMatrix,I,J} |
| 17 | +# prob::OP |
| 18 | +# √Σ::S0 |
| 19 | +# Xp::X |
| 20 | +# u_idx::I |
| 21 | +# θ_idx::J |
| 22 | +# function UncertainODEProblem(prob::OP) where OP<:ODEProblem |
| 23 | +# u_idx = findall(isuncertain, prob.u0) |
| 24 | +# θ_idx = findall(isuncertain, prob.θ) |
| 25 | +# |
| 26 | +# uσ = [u.σ for u in prob.u0[u_idx]] |
| 27 | +# θσ = [θ.σ for θ in prob.θ[θ_idx]] |
| 28 | +# √Σ = Diagonal(vcat(uσ,θσ)) |
| 29 | +# |
| 30 | +# n = (length(uσ)+length(θσ))*2 |
| 31 | +# Qp = sqrt(n) * [I(n) -I(n)] |
| 32 | +# Xp = vcat(prob.u0,prob.θ) .+ √Σ*Qp |
| 33 | +# Σ = √Σ * √Σ' |
| 34 | +# |
| 35 | +# new{OP,typeof(√Σ),typeof(u_idx),typeof(θ_idx)}(prob,√Σ,u_idx,θ_idx) |
| 36 | +# end |
| 37 | +#end |
| 38 | +struct UncertainODEProblem{OP,US,I,J} <: AbstractODEProblem |
| 39 | + prob::OP |
| 40 | + U0::US |
| 41 | + u_idx::I |
| 42 | + θ_idx::J |
| 43 | + function UncertainODEProblem(prob::ODEProblem) |
| 44 | + u_idx = findall(isuncertain, prob.u0) |
| 45 | + θ_idx = findall(isuncertain, prob.θ) |
| 46 | + idx = vcat(u_idx, length(prob.u0) .+ θ_idx) |
| 47 | + |
| 48 | + n = length(idx) |
| 49 | + Qp = hcat(map(idx) do i |
| 50 | + q = zeros(length(prob.u0)+length(prob.θ)) |
| 51 | + q[i] = sqrt(n) |
| 52 | + q |
| 53 | + end, map(idx) do i |
| 54 | + q = zeros(length(prob.u0)+length(prob.θ)) |
| 55 | + q[i] = -sqrt(n) |
| 56 | + q |
| 57 | + end) |
| 58 | + Qp = reduce(hcat,Qp) |
| 59 | + Σ_ = Diagonal(vcat(sig.(prob.u0), sig.(prob.θ))) |
| 60 | + μ0 = vcat(mu.(prob.u0),mu.(prob.θ)) .+ Σ_*Qp |
| 61 | + #U0 = μ0 .± diag(Σ_) |
| 62 | + |
| 63 | + prob = ODEProblem(prob.f, prob.tspan, mu.(prob.u0), mu.(prob.θ)) |
| 64 | + |
| 65 | + new{typeof(prob),typeof(μ0),typeof(u_idx),typeof(θ_idx)}(prob,μ0,u_idx,θ_idx) |
| 66 | + end |
| 67 | +end |
| 68 | + |
| 69 | +nr_uncertainties(p::UncertainODEProblem) = length(p.u_idx)+length(p.θ_idx) |
| 70 | + |
| 71 | +struct UncertainODESolver{S<:ODESolver} <: ODESolver |
| 72 | + solver::S |
| 73 | +end |
| 74 | + |
| 75 | +# function get_xp(p::UncertainODEProblem,i::Int) |
| 76 | +# xp = p.Xp[:,i] |
| 77 | +# u = xp[p.u_idx] |
| 78 | +# θ = p.prob.θ |
| 79 | +# θ[p.θ_idx] .= xp[length(p.u_idx)+1:end] |
| 80 | +# (u,θ) |
| 81 | +# end |
| 82 | + |
| 83 | + |
| 84 | +function setmean!(p::UncertainODEProblem, x) |
| 85 | + nu = length(p.prob.u0) |
| 86 | + p.prob.θ .= x[nu .+ (1:length(p.prob.θ))] |
| 87 | + u = x[1:nu] |
| 88 | +end |
| 89 | + |
| 90 | +function (s::UncertainODESolver)(p::UncertainODEProblem, μs, t) |
| 91 | + N = nr_uncertainties(p) |
| 92 | + μs = map(1:N) do i |
| 93 | + u = setmean!(p, μs[:,i]) |
| 94 | + u = s.solver(p.prob, u, t)[1] |
| 95 | + vcat(u, p.prob.θ) |
| 96 | + end |
| 97 | + μs = reduce(hcat, μs) |
| 98 | + μ = mean(μs,dims=2) |
| 99 | + Σ = Matrix((μs .- μ)*(μs .- μ)'/N) |
| 100 | + σ = sqrt.(diag(Σ)) |
| 101 | + σ[p.u_idx] .= 0 |
| 102 | + σ[p.θ_idx .+ length(p.prob.u0)] .= 0 |
| 103 | + #μ .± σ, t+1 |
| 104 | + μs, t+1 |
| 105 | +end |
| 106 | + |
| 107 | +function solve(p::UncertainODEProblem, solver::UncertainODESolver) |
| 108 | + t = p.prob.tspan[1]; u = p.U0 |
| 109 | + us = [u]; ts = [t] |
| 110 | + while t < prob.tspan[2] |
| 111 | + (u,t) = solver(p, u, t) |
| 112 | + push!(us,u) |
| 113 | + push!(ts,t) |
| 114 | + end |
| 115 | + ts, reduce(hcat,us) |
| 116 | +end |
| 117 | + |
| 118 | + |
| 119 | + |
| 120 | +uprob = UncertainODEProblem(prob) |
| 121 | +solver = UncertainODESolver(RK2(0.2)) |
| 122 | +t, X = solve(uprob,solver) |
| 123 | + |
| 124 | +# function solve(f,x0::AbstractVector,sqΣ0, θ,dt,N) |
| 125 | +# n = length(x0) |
| 126 | +# n2 = 2*length(x0) |
| 127 | +# Qp = sqrt(n)*[I(n) -I(n)] |
| 128 | +# |
| 129 | +# X = hcat([zero(x0) for i=1:N]...) |
| 130 | +# S = hcat([zero(x0) for i=1:N]...) |
| 131 | +# X[:,1]=x0 |
| 132 | +# Xp = x0 .+ sqΣ0*Qp |
| 133 | +# sqΣ = sqΣ0 |
| 134 | +# Σ = sqΣ* sqΣ' |
| 135 | +# S[:,1]= diag(Σ) |
| 136 | +# for t=1:N-1 |
| 137 | +# for i=1:n2 # all quadrature points |
| 138 | +# Xp[:,i].=Xp[:,i] + dt*f(Xp[:,i],θ) |
| 139 | +# end |
| 140 | +# mXp=mean(Xp,dims=2) |
| 141 | +# X[:,t+1]=mXp |
| 142 | +# Σ=Matrix((Xp.-mXp)*(Xp.-mXp)'/n2) |
| 143 | +# S[:,t+1]=sqrt.(diag(Σ)) |
| 144 | +# # @show Σ |
| 145 | +# |
| 146 | +# end |
| 147 | +# X,S |
| 148 | +# end |
| 149 | + |
| 150 | + |
| 151 | + |
| 152 | + |
| 153 | + |
| 154 | + |
| 155 | + |
| 156 | + |
0 commit comments