sdptest.jl

using JuMP, Mosek, FactCheck

sdp_solvers = [MosekSolver(LOG=0)]
include(Pkg.dir("JuMP","test","data","robust_uncertainty.jl"))

facts("[sdp] Robust uncertainty example") do
for solver in sdp_solvers
context("With solver $(typeof(solver))") do
    R = 1
    d = 3
    𝛿 = 0.05
    ɛ = 0.05
    N = ceil((2+2log(2/𝛿))^2) + 1

    Γ1(𝛿,N) = (R/sqrt(N))*(2+sqrt(2*log(1/𝛿)))
    Γ2(𝛿,N) = (2R^2/sqrt(N))*(2+sqrt(2*log(2/𝛿)))

    for d in [3]#,5,8]; context("d = $d") do

        μhat = μhats[d]
        M = Ms[d]
        Σhat = 1/(d-1)*(M-ones(d)*μhat')'*(M-ones(d)*μhat')

        m = Model(solver=solver)

        @defVar(m, Σ[1:d,1:d], SDP)
        @defVar(m, u[1:d])
        @defVar(m, μ[1:d])

        @defVar(m, t1 >= 0)
        @defVar(m, L1[1:d])
        @addConstraint(m, L1 .== (μ-μhat))
        @addConstraint(m, sum{L1[i]^2, i=1:d} <= t1^2)
        @addConstraint(m, t1 <= Γ1(𝛿/2,N))

        @defVar(m, t2 >= 0)
        @defVar(m, L2[1:d,1:d])
        @addConstraint(m, L2 .== (Σ-Σhat))
        @addConstraint(m, sum{L2[i,j]^2, i=1:d, j=1:d} <= t2^2)
        @addConstraint(m, t2 <= Γ2(𝛿/2,N))

        A = [(1-ɛ)/ɛ (u-μ)';
             (u-μ)     Σ   ]
        @addSDPConstraint(m, A >= 0)

        c = cs[d]
        @setObjective(m, Max, dot(c,u))

        stat = solve(m)

        object = getObjectiveValue(m)
        exact = dot(μhat,c) + Γ1(𝛿/2,N)*norm(c) + sqrt((1-ɛ)/ɛ)*sqrt(dot(c,(Σhat+Γ2(𝛿/2,N)*eye(d,d))*c))
        @show exact
        @show object
        @fact stat --> :Optimal
        @fact abs(object - exact) --> roughly(0, 1e-5)
    end; #end
end; end; end