Registering the user-defined function "multivariate normal cdf" for nonlinear modeling in JuMP

Specific Domains Optimization (Mathematical)
,
3

You can set up prob(x,y) as the first value returned from mvnormcdf:

prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)[1]

Also, your cdf in ∇prob need the first argument to work, e.g.:

cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2))) #I added Normal(), not sure if this is correct for you

Then this works:


using JuMP, LinearAlgebra
using Distributions, NLopt
using Random
using MvNormalCDF

begin
    ex2 = Model(NLopt.Optimizer)
    set_optimizer_attribute(ex2, "algorithm", :LD_MMA)

    @variable(ex2, x[1:2] >= 0)

    @objective(ex2, Min, x[1] + x[2])

    set_start_value(x[1], 5)
    set_start_value(x[2], 5)

    p = 0.9 # probability

    μ = [5; 5] # mean
    vari = [4; 4] # variance
    ρ = 0.5 # correlation
    Σ = [vari[1] ρ; ρ vari[2]] # covariance matrix
    a = [-Inf; -Inf] 

    norm_cdf(x) = cdf(Normal(), x)
    norm_pdf(x) = exp(-(x^2) / 2) / sqrt(2 * pi)
    register(ex2, :cdf, 1, norm_cdf, norm_pdf; autodiff = true)

    # Define the mvnormcdf and its gradient, then register the function 
    prob(x, y) = mvnormcdf(Σ, a-μ, [x; y]-μ)[1]
    function ∇prob(g::AbstractVector{T}, x::T, y::T) where {T}
        g[1] = norm_pdf((x - μ[1])/sqrt(vari[1])) * cdf(Normal(),(y-(μ[2]+ρ/sqrt(vari[1])*(x-μ[1])))/sqrt(vari[2]^2*(1- ρ^2)))
        g[2] = norm_pdf((y - μ[2])/sqrt(vari[2]))* cdf(Normal(),(x-(μ[1]+ρ/sqrt(vari[2])*(y-μ[2])))/sqrt(vari[1]^2*(1- ρ^2)))
        return
    end


    register(ex2, :mvncdf, 2, prob, ∇prob)


    # Use the mvncdf in the nonlinear constraint
    # This is where the type error occurs
    @NLconstraint(ex2,mvncdf(x[1], x[2]) >= p)

    JuMP.optimize!(ex2)

    println(value.(x))

end

And I got:

[8.247724657052249, 8.247724657052249]
3 Likes