How to best handle these exceptions in FOC of maximization?

I’m trying to maximize a function f(t_1, t_2) for which the first order conditions yield 2 equations depending on some conditions. The FOCs divide the plane into 3 regions, where t_1 > t_2, where t_1 < t_2, and where t_1 = t_2. I am solving each of those systems of equations for each of those cases, but there are cases that are not feasible for given parameters, so I have to handle possible errors when I feed the solver a system that is not feasible.

I am struggling to find a MWE, but I hope this example can convey the essence of the problem. For simplicity, let us consider only the regions where t_1 > t_2, and where t_1 < t_2.

using Optim
using Roots

function H(s; θ=θ)
    s^2 - θ

function N(t1, t2; θ=θ)
    -((t1 - 10)^2 + (t2 - 10)^2) + θ

function foc_t1t2(t2bound; θ=θ)
    t1star = find_zero(t1 -> H(t1; θ=θ), (0.0, 100_000.0))

    solt2 = optimize(t2 -> -N(t1star, t2; θ=θ), t1star, t2bound, Optim.Brent())
    t2star = solt2.minimizer

    return [t1star, t2star]

function foc_t2t1(t1bound; θ=θ)
    t2star = find_zero(t2 -> H(t2; θ=θ), (0.0, 100_000.0))

    solt1 = optimize(t1 -> -N(t1, t2star; θ=θ), t2star, t1bound, Optim.Brent())
    t1star = solt1.minimizer

    return [t1star, t2star]

function sol_foc(; θ1=θ1, θ2=θ2, K1=K1, K2=K2)

    # bounds
    H1inv(y1) = Roots.find_zero(t1 -> H(t1; θ=θ1) + K1 - y1, (0.0, 100_000.0))
    t1bound = H1inv(K1)
    H2inv(y2) = Roots.find_zero(t2 -> H(t2; θ=θ2) + K2 - y2, (0.0, 100_000.0))
    t2bound = H2inv(K2)

    sol_t1t2 = [NaN, NaN]
        sol_t1t2 .= foc_t1t2(t2bound; θ=θ1)
        @info "t1 < t2 not feasible"

    sol_t2t1 = [NaN, NaN]
        sol_t2t1 .= foc_t2t1(t1bound; θ=θ2)
        @info "t2 < t1 not feasible"

    # test which one yields largest value
    trueval = -Inf
    truezero = [NaN, NaN]
    for v in (sol_t1t2, sol_t2t1)
        n = N(v...; θ=θ1+θ2)
        if trueval < n
            trueval = n
            truezero .= v

    return truezero

θ1 = 0
θ2 = 2
K1 = 0
K2 = 100

those parameters then yield this

julia> sol_foc(; θ1=θ1, θ2=θ2, K1=K1, K2=K2)
[ Info: t2 < t1 not feasible
2-element Vector{Float64}:

I’m looking to see if there is a better (more performant) way of dealing with those exceptions.

Hey @amrods,

without understanding anything about your problem domain I checked the MWE with @btime and it showed

141.600 μs (122 allocations: 5.33 KiB)

so my only proposal would be to run them in parallel (not much of an improvement and only working if your tasks are big enough, I know)?

1 Like