|
| 1 | +@kwdef struct AugLag |
| 2 | + inner::Any |
| 3 | + τ = 0.5 |
| 4 | + γ = 10.0 |
| 5 | + λmin = -1e20 |
| 6 | + λmax = 1e20 |
| 7 | + μmin = 0.0 |
| 8 | + μmax = 1e20 |
| 9 | + ϵ = 1e-8 |
| 10 | +end |
| 11 | + |
| 12 | +SciMLBase.supports_opt_cache_interface(::AugLag) = true |
| 13 | +SciMLBase.allowsbounds(::AugLag) = true |
| 14 | +SciMLBase.requiresgradient(::AugLag) = true |
| 15 | +SciMLBase.allowsconstraints(::AugLag) = true |
| 16 | +SciMLBase.requiresconsjac(::AugLag) = true |
| 17 | + |
| 18 | +function __map_optimizer_args(cache::Optimization.OptimizationCache, opt::AugLag; |
| 19 | + callback = nothing, |
| 20 | + maxiters::Union{Number, Nothing} = nothing, |
| 21 | + maxtime::Union{Number, Nothing} = nothing, |
| 22 | + abstol::Union{Number, Nothing} = nothing, |
| 23 | + reltol::Union{Number, Nothing} = nothing, |
| 24 | + verbose::Bool = false, |
| 25 | + kwargs...) |
| 26 | + if !isnothing(abstol) |
| 27 | + @warn "common abstol is currently not used by $(opt)" |
| 28 | + end |
| 29 | + if !isnothing(maxtime) |
| 30 | + @warn "common abstol is currently not used by $(opt)" |
| 31 | + end |
| 32 | + |
| 33 | + mapped_args = (;) |
| 34 | + |
| 35 | + if cache.lb !== nothing && cache.ub !== nothing |
| 36 | + mapped_args = (; mapped_args..., lb = cache.lb, ub = cache.ub) |
| 37 | + end |
| 38 | + |
| 39 | + if !isnothing(maxiters) |
| 40 | + mapped_args = (; mapped_args..., maxiter = maxiters) |
| 41 | + end |
| 42 | + |
| 43 | + if !isnothing(reltol) |
| 44 | + mapped_args = (; mapped_args..., pgtol = reltol) |
| 45 | + end |
| 46 | + |
| 47 | + return mapped_args |
| 48 | +end |
| 49 | + |
| 50 | +function SciMLBase.__solve(cache::OptimizationCache{ |
| 51 | + F, |
| 52 | + RC, |
| 53 | + LB, |
| 54 | + UB, |
| 55 | + LC, |
| 56 | + UC, |
| 57 | + S, |
| 58 | + O, |
| 59 | + D, |
| 60 | + P, |
| 61 | + C |
| 62 | +}) where { |
| 63 | + F, |
| 64 | + RC, |
| 65 | + LB, |
| 66 | + UB, |
| 67 | + LC, |
| 68 | + UC, |
| 69 | + S, |
| 70 | + O <: |
| 71 | + AugLag, |
| 72 | + D, |
| 73 | + P, |
| 74 | + C |
| 75 | +} |
| 76 | + maxiters = Optimization._check_and_convert_maxiters(cache.solver_args.maxiters) |
| 77 | + |
| 78 | + local x |
| 79 | + |
| 80 | + solver_kwargs = __map_optimizer_args(cache, cache.opt; maxiters, cache.solver_args...) |
| 81 | + |
| 82 | + if !isnothing(cache.f.cons) |
| 83 | + eq_inds = [cache.lcons[i] == cache.ucons[i] for i in eachindex(cache.lcons)] |
| 84 | + ineq_inds = (!).(eq_inds) |
| 85 | + |
| 86 | + τ = cache.opt.τ |
| 87 | + γ = cache.opt.γ |
| 88 | + λmin = cache.opt.λmin |
| 89 | + λmax = cache.opt.λmax |
| 90 | + μmin = cache.opt.μmin |
| 91 | + μmax = cache.opt.μmax |
| 92 | + ϵ = cache.opt.ϵ |
| 93 | + |
| 94 | + λ = zeros(eltype(cache.u0), sum(eq_inds)) |
| 95 | + μ = zeros(eltype(cache.u0), sum(ineq_inds)) |
| 96 | + |
| 97 | + cons_tmp = zeros(eltype(cache.u0), length(cache.lcons)) |
| 98 | + cache.f.cons(cons_tmp, cache.u0) |
| 99 | + ρ = max(1e-6, |
| 100 | + min(10, 2 * (abs(cache.f(cache.u0, iterate(cache.p)[1]))) / norm(cons_tmp))) |
| 101 | + |
| 102 | + _loss = function (θ, p = cache.p) |
| 103 | + x = cache.f(θ, p) |
| 104 | + cons_tmp .= zero(eltype(θ)) |
| 105 | + cache.f.cons(cons_tmp, θ) |
| 106 | + cons_tmp[eq_inds] .= cons_tmp[eq_inds] - cache.lcons[eq_inds] |
| 107 | + cons_tmp[ineq_inds] .= cons_tmp[ineq_inds] .- cache.ucons[ineq_inds] |
| 108 | + opt_state = Optimization.OptimizationState(u = θ, objective = x[1]) |
| 109 | + if cache.callback(opt_state, x...) |
| 110 | + error("Optimization halted by callback.") |
| 111 | + end |
| 112 | + return x[1] + sum(@. λ * cons_tmp[eq_inds] + ρ / 2 * (cons_tmp[eq_inds] .^ 2)) + |
| 113 | + 1 / (2 * ρ) * sum((max.(Ref(0.0), μ .+ (ρ .* cons_tmp[ineq_inds]))) .^ 2) |
| 114 | + end |
| 115 | + |
| 116 | + prev_eqcons = zero(λ) |
| 117 | + θ = cache.u0 |
| 118 | + β = max.(cons_tmp[ineq_inds], Ref(0.0)) |
| 119 | + prevβ = zero(β) |
| 120 | + eqidxs = [eq_inds[i] > 0 ? i : nothing for i in eachindex(ineq_inds)] |
| 121 | + ineqidxs = [ineq_inds[i] > 0 ? i : nothing for i in eachindex(ineq_inds)] |
| 122 | + eqidxs = eqidxs[eqidxs .!= nothing] |
| 123 | + ineqidxs = ineqidxs[ineqidxs .!= nothing] |
| 124 | + function aug_grad(G, θ, p) |
| 125 | + cache.f.grad(G, θ, p) |
| 126 | + if !isnothing(cache.f.cons_jac_prototype) |
| 127 | + J = Float64.(cache.f.cons_jac_prototype) |
| 128 | + else |
| 129 | + J = zeros((length(cache.lcons), length(θ))) |
| 130 | + end |
| 131 | + cache.f.cons_j(J, θ) |
| 132 | + __tmp = zero(cons_tmp) |
| 133 | + cache.f.cons(__tmp, θ) |
| 134 | + __tmp[eq_inds] .= __tmp[eq_inds] .- cache.lcons[eq_inds] |
| 135 | + __tmp[ineq_inds] .= __tmp[ineq_inds] .- cache.ucons[ineq_inds] |
| 136 | + G .+= sum( |
| 137 | + λ[i] .* J[idx, :] + ρ * (__tmp[idx] .* J[idx, :]) |
| 138 | + for (i, idx) in enumerate(eqidxs); |
| 139 | + init = zero(G)) #should be jvp |
| 140 | + G .+= sum( |
| 141 | + 1 / ρ * (max.(Ref(0.0), μ[i] .+ (ρ .* __tmp[idx])) .* J[idx, :]) |
| 142 | + for (i, idx) in enumerate(ineqidxs); |
| 143 | + init = zero(G)) #should be jvp |
| 144 | + end |
| 145 | + |
| 146 | + opt_ret = ReturnCode.MaxIters |
| 147 | + n = length(cache.u0) |
| 148 | + |
| 149 | + augprob = OptimizationProblem( |
| 150 | + OptimizationFunction(_loss; grad = aug_grad), cache.u0, cache.p) |
| 151 | + |
| 152 | + solver_kwargs = Base.structdiff(solver_kwargs, (; lb = nothing, ub = nothing)) |
| 153 | + |
| 154 | + for i in 1:(maxiters / 10) |
| 155 | + prev_eqcons .= cons_tmp[eq_inds] .- cache.lcons[eq_inds] |
| 156 | + prevβ .= copy(β) |
| 157 | + res = solve(augprob, cache.opt.inner, maxiters = maxiters / 10) |
| 158 | + θ = res.u |
| 159 | + cons_tmp .= 0.0 |
| 160 | + cache.f.cons(cons_tmp, θ) |
| 161 | + λ = max.(min.(λmax, λ .+ ρ * (cons_tmp[eq_inds] .- cache.lcons[eq_inds])), λmin) |
| 162 | + β = max.(cons_tmp[ineq_inds], -1 .* μ ./ ρ) |
| 163 | + μ = min.(μmax, max.(μ .+ ρ * cons_tmp[ineq_inds], μmin)) |
| 164 | + if max(norm(cons_tmp[eq_inds] .- cache.lcons[eq_inds], Inf), norm(β, Inf)) > |
| 165 | + τ * max(norm(prev_eqcons, Inf), norm(prevβ, Inf)) |
| 166 | + ρ = γ * ρ |
| 167 | + end |
| 168 | + if norm( |
| 169 | + (cons_tmp[eq_inds] .- cache.lcons[eq_inds]) ./ cons_tmp[eq_inds], Inf) < |
| 170 | + ϵ && norm(β, Inf) < ϵ |
| 171 | + opt_ret = ReturnCode.Success |
| 172 | + break |
| 173 | + end |
| 174 | + end |
| 175 | + stats = Optimization.OptimizationStats(; iterations = maxiters, |
| 176 | + time = 0.0, fevals = maxiters, gevals = maxiters) |
| 177 | + return SciMLBase.build_solution( |
| 178 | + cache, cache.opt, θ, x, |
| 179 | + stats = stats, retcode = opt_ret) |
| 180 | + end |
| 181 | +end |
0 commit comments