Skip to content

Commit

Permalink
don't automatically transform quadratic constraints into conic form
Browse files Browse the repository at this point in the history
  • Loading branch information
@mlubin
mlubin committed Mar 25, 2017
1 parent 1f0fbaf commit 3e66b34
Show file tree
Hide file tree
Showing 4 changed files with 72 additions and 166 deletions.
123 changes: 4 additions & 119 deletions src/solvers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -282,60 +282,6 @@ function solve(m::Model; suppress_warnings=false,
stat
end

function isquadsoc(m::Model)
# check if all quadratic constraints are actually conic
all_conic = true
for qconstr in m.quadconstr
q = copy(qconstr.terms)
if qconstr.sense == :(>=)
q *= -1
end
if !(all(t->t==0, q.aff.coeffs) && q.aff.constant == 0)
all_conic = false
break
end
n_pos_on_diag = 0
off_diag_idx = 0
neg_diag_idx = 0
n = length(q.qvars1)
nz = 0
for i in 1:n
q.qcoeffs[i] == 0 && continue
nz += 1
if q.qvars1[i].col == q.qvars2[i].col
if q.qcoeffs[i] == 1
n_pos_on_diag += 1
elseif q.qcoeffs[i] == -1
if !(neg_diag_idx == off_diag_idx == 0)
all_conic = false; break
end
neg_diag_idx = i
else
all_conic = false; break
end
else
if q.qcoeffs[i] == -1
if !(neg_diag_idx == off_diag_idx == 0)
all_conic = false; break
end
off_diag_idx = i
else
all_conic = false; break
end
end
end
if n_pos_on_diag == nz-1 && neg_diag_idx > 0
# Plain SOC
elseif n_pos_on_diag == nz-1 && off_diag_idx > 0
# Rotated SOC
else
all_conic = false
end
end
return all_conic && length(m.quadconstr) > 0

end

# Converts the JuMP Model into a MathProgBase model based on the
# traits of the model
function build(m::Model; suppress_warnings=false, relaxation=false, traits=ProblemTraits(m,relaxation=relaxation))
Expand All @@ -347,11 +293,6 @@ function build(m::Model; suppress_warnings=false, relaxation=false, traits=Probl
# to build the problem
traits.nlp && return _buildInternalModel_nlp(m, traits)

## Temporary hack for Mosek, see https://github.com/JuliaOpt/Mosek.jl/pull/67
if contains("$(typeof(m.solver))","MosekSolver") && isquadsoc(m)
traits.conic = true
end

if traits.conic
# If there are semicontinuous/semi-integer variables, we will have to
# adjust the b vector below to construct a valid relaxation. This seems
Expand All @@ -361,6 +302,7 @@ function build(m::Model; suppress_warnings=false, relaxation=false, traits=Probl
end

traits.qp && error("JuMP does not support quadratic objectives for conic problems")
traits.qc && error("JuMP does not support mixing quadratic and conic constraints")

# Obtain a fresh MPB model for the solver
# If the problem is conic, we rebuild the problem from
Expand Down Expand Up @@ -696,63 +638,6 @@ function conicdata(m::Model)

soc_cones = Any[]
rsoc_cones = Any[]
numQuadRows = 0
for qconstr in m.quadconstr
q = copy(qconstr.terms)
if qconstr.sense == :(>=)
q *= -1
end
if !(all(t->t==0, q.aff.coeffs) && q.aff.constant == 0)
error("Quadratic constraint $qconstr must be in second-order cone form")
end
n_pos_on_diag = 0
off_diag_idx = 0
neg_diag_idx = 0
n = length(q.qvars1)
nz = 0
for i in 1:n
q.qcoeffs[i] == 0 && continue
nz += 1
if q.qvars1[i].col == q.qvars2[i].col
if q.qcoeffs[i] == 1
n_pos_on_diag += 1
elseif q.qcoeffs[i] == -1
neg_diag_idx == off_diag_idx == 0 || error("Invalid SOC constraint $qconstr")
neg_diag_idx = i
end
else
if q.qcoeffs[i] == -1
neg_diag_idx == off_diag_idx == 0 || error("Invalid rotated SOC constraint $qconstr")
off_diag_idx = i
nz += 1
end
end
end
cone = Array{Int}(nz)
if n_pos_on_diag == nz-1 && neg_diag_idx > 0
cone[1] = q.qvars1[neg_diag_idx].col
r = 1
for i in 1:n
(q.qcoeffs[i] == 0 || i == neg_diag_idx) && continue
r += 1
cone[r] = q.qvars1[i].col
end
push!(soc_cones, cone)
elseif n_pos_on_diag == nz-2 && off_diag_idx > 0
cone[1] = q.qvars1[off_diag_idx].col
cone[2] = q.qvars2[off_diag_idx].col
r = 2
for i in 1:n
(q.qcoeffs[i] == 0 || i == off_diag_idx) && continue
r += 1
cone[r] = q.qvars1[i].col
end
push!(rsoc_cones, cone)
else
error("Quadratic constraint $qconstr is not conic representable")
end
numQuadRows += length(cone)
end

linconstr = m.linconstr::Vector{LinearConstraint}
numLinRows = length(linconstr)
Expand Down Expand Up @@ -815,7 +700,7 @@ function conicdata(m::Model)
numNormRows += 1
numSOCRows += length(con.normexpr.norm.terms) + 1
end
numRows = numLinRows + numBounds + numQuadRows + numSOCRows + numSDPRows + numSymRows
numRows = numLinRows + numBounds + numSOCRows + numSDPRows + numSymRows

# should maintain the order of constraints in the above form
# throughout the code c is the conic constraint index
Expand Down Expand Up @@ -931,7 +816,7 @@ function conicdata(m::Model)
b[rng] = 0
c += n
end
@assert c == numLinRows + numBounds + numQuadRows
@assert c == numLinRows + numBounds

tmpelts = tmprow.elts
tmpnzidx = tmprow.nzidx
Expand Down Expand Up @@ -961,7 +846,7 @@ function conicdata(m::Model)
push!(con_cones, (:SOC, soc_start:c))
constr_dual_map[numLinRows + numBounds + socidx] = collect(soc_start:c)
end
@assert c == numLinRows + numBounds + numQuadRows + numSOCRows
@assert c == numLinRows + numBounds + numSOCRows

sdpconstr_sym = Vector{Vector{Tuple{Int,Int}}}(length(m.sdpconstr))
numDroppedSym = 0
Expand Down
73 changes: 38 additions & 35 deletions test/qcqpmodel.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,6 @@ using Base.Test
@test MathProgBase.numquadconstr(modQ) == 1
@test MathProgBase.numlinconstr(modQ) == 0
@test MathProgBase.numconstr(modQ) == 1
@test JuMP.isquadsoc(modQ) == false

@test solve(modQ) == :Optimal
@test isapprox(getobjectivevalue(modQ), -1-4/sqrt(3), atol=1e-6)
Expand Down Expand Up @@ -186,14 +185,47 @@ using Base.Test

@testset "SOC constraints (continuous) with $solver" for solver in soc_solvers

# SOC version 1
modN = Model(solver=solver)
@variable(modN, x)
@variable(modN, y)
@variable(modN, t >= 0)
@objective(modN, Min, t)
@constraint(modN, x + y >= 1)
@constraint(modN, norm([x,y]) <= t)

# Getter/setters
@test JuMP.numsocconstr(modN) == 1
@test MathProgBase.numconstr(modN) == 2

@test solve(modN) == :Optimal
@test isapprox(getobjectivevalue(modN), sqrt(1/2), atol=1e-6)
@test isapprox([getvalue(x), getvalue(y), getvalue(t)], [0.5,0.5,sqrt(1/2)], atol=1e-3)

# SOC version 2
modN = Model(solver=solver)
@variable(modN, x)
@variable(modN, y)
@variable(modN, t >= 0)
@objective(modN, Min, t)
@constraint(modN, x + y >= 1)
tmp = [x,y]
@constraint(modN, norm(tmp[i] for i=1:2) <= t)

@test solve(modN) == :Optimal
@test isapprox(getobjectivevalue(modN), sqrt(1/2), atol=1e-6)
@test isapprox([getvalue(x), getvalue(y), getvalue(t)], [0.5,0.5,sqrt(1/2)], atol=1e-3)
end

@testset "SOC constraints (continuous) in quadratic form with $solver" for solver in quad_soc_solvers

modQ = Model(solver=solver)
@variable(modQ, x)
@variable(modQ, y)
@variable(modQ, t >= 0)
@objective(modQ, Min, t)
@constraint(modQ, x+y >= 1)
@constraint(modQ, x^2 + y^2 <= t^2)
@test JuMP.isquadsoc(modQ) == true

@test solve(modQ) == :Optimal
@test isapprox(getobjectivevalue(modQ), sqrt(1/2), atol=1e-6)
Expand Down Expand Up @@ -230,37 +262,6 @@ using Base.Test
@test solve(modV) == :Optimal
@test isapprox(getobjectivevalue(modV), sqrt(1/2), atol=1e-6)
@test isapprox([getvalue(x), getvalue(y), getvalue(t)], [0.5,0.5,sqrt(1/2)], atol=1e-3)

# SOC version 1
modN = Model(solver=solver)
@variable(modN, x)
@variable(modN, y)
@variable(modN, t >= 0)
@objective(modN, Min, t)
@constraint(modN, x + y >= 1)
@constraint(modN, norm([x,y]) <= t)

# Getter/setters
@test JuMP.numsocconstr(modN) == 1
@test MathProgBase.numconstr(modN) == 2

@test solve(modN) == :Optimal
@test isapprox(getobjectivevalue(modN), sqrt(1/2), atol=1e-6)
@test isapprox([getvalue(x), getvalue(y), getvalue(t)], [0.5,0.5,sqrt(1/2)], atol=1e-3)

# SOC version 2
modN = Model(solver=solver)
@variable(modN, x)
@variable(modN, y)
@variable(modN, t >= 0)
@objective(modN, Min, t)
@constraint(modN, x + y >= 1)
tmp = [x,y]
@constraint(modN, norm(tmp[i] for i=1:2) <= t)

@test solve(modN) == :Optimal
@test isapprox(getobjectivevalue(modN), sqrt(1/2), atol=1e-6)
@test isapprox([getvalue(x), getvalue(y), getvalue(t)], [0.5,0.5,sqrt(1/2)], atol=1e-3)
end

@testset "SOC duals with $solver" for solver in soc_solvers
Expand Down Expand Up @@ -365,16 +366,18 @@ using Base.Test
@test isapprox(getobjectivevalue(modQ), -0.75, atol=1e-6)
end

@testset "Rotated second-order cones with $solver" for solver in rsoc_solvers
@testset "Rotated second-order cones with $solver" for solver in quad_soc_solvers
mod = Model(solver=solver)

@variable(mod, x[1:5] >= 0)
@variable(mod, 0 <= u <= 5)
@variable(mod, v)
@variable(mod, t1 == 1)
@variable(mod, t2 == 1)

@objective(mod, Max, v)

@constraint(mod, norm(x) <= 1)
@constraint(mod, sum(x.^2) <= t1*t2)
@constraint(mod, v^2 <= u * x[1])

@test solve(mod) == :Optimal
Expand Down
37 changes: 25 additions & 12 deletions test/sdp.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,12 +150,11 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
@test macroexpand(:(@variable(m, -rand(5,5) <= nonsymmetric[1:5,1:5] <= rand(5,5), Symmetric))).head == :error
end

@testset "SDP with quadratics with $solver" for solver in sdp_solvers
@testset "SDP with SOC with $solver" for solver in sdp_solvers
m = Model(solver=solver)
@variable(m, X[1:2,1:2], SDP)
@variable(m, y[0:2])
@constraint(m, y[0] >= 0)
@constraint(m, y[1]^2 + y[2]^2 <= y[0]^2)
@constraint(m, norm([y[1],y[2]]) <= y[0])
@SDconstraint(m, X <= eye(2))
@constraint(m, X[1,1] + X[1,2] == y[1] + y[2])
@objective(m, Max, trace(X) - y[0])
Expand Down Expand Up @@ -463,17 +462,13 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
@variable(m, u[1:d])
@variable(m, μ[1:d])

@variable(m, t1 >= 0)
@variable(m, L1[1:d])
@constraint(m, L1 .== (μ-μhat))
@constraint(m, sum(L1[i]^2 for i=1:d) <= t1^2)
@constraint(m, t1 <= Γ1(𝛿/2,N))
@constraint(m, norm(L1) <= Γ1(𝛿/2,N))

@variable(m, t2 >= 0)
@variable(m, L2[1:d,1:d])
@constraint(m, L2 .== (Σ-Σhat))
@constraint(m, sum(L2[i,j]^2 for i=1:d, j=1:d) <= t2^2)
@constraint(m, t2 <= Γ2(𝛿/2,N))
@constraint(m, vecnorm(L2) <= Γ2(𝛿/2,N))

A = [(1-ɛ)/ɛ (u-μ)';
(u-μ) Σ ]
Expand Down Expand Up @@ -636,7 +631,13 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
@test all(isnan(getdual(c)))
status = solve(m)

@test status == :Optimal
if contains(string(typeof(solver)),"MosekSolver")
# Mosek returns Stall on this instance
# Hack until we fix statuses in MPB
JuMP.fillConicDuals(m)
else
@test status == :Optimal
end
@test isapprox(getobjectivevalue(m), 0, atol=1e-5)
@test isapprox(getvalue(y), 0, atol=1e-5)

Expand Down Expand Up @@ -680,7 +681,13 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
@objective(m, Max, y/2+z/2)
status = solve(m)

@test status == :Optimal
if contains(string(typeof(solver)),"MosekSolver")
# Mosek returns Stall on this instance
# Hack until we fix statuses in MPB
JuMP.fillConicDuals(m)
else
@test status == :Optimal
end
@test isapprox(getobjectivevalue(m), 0, atol=1e-5)
@test isapprox(getvalue(y), 0, atol=1e-5)
@test isapprox(getvalue(z), 0, atol=1e-5)
Expand All @@ -706,7 +713,13 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
@objective(m, Max, y)
status = solve(m)

@test status == :Optimal
if contains(string(typeof(solver)),"MosekSolver")
# Mosek returns Stall on this instance
# Hack until we fix statuses in MPB
JuMP.fillConicDuals(m)
else
@test status == :Optimal
end
@test isapprox(getobjectivevalue(m), 0, atol=1e-5)
@test isapprox(getvalue(y), 0, atol=1e-5)
@test isapprox(getvalue(z), 0, atol=1e-5)
Expand Down
5 changes: 5 additions & 0 deletions test/solvers.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,11 @@ cpx && push!(quad_solvers, CPLEX.CplexSolver(CPX_PARAM_SCRIND=0))
xpr && push!(quad_solvers, Xpress.XpressSolver(OUTPUTLOG=0, FEASTOL = 1e-9, BARPRIMALSTOP = 1e-9, BARGAPSTOP = 1e-9, BARDUALSTOP = 1e-9))
mos && push!(quad_solvers, Mosek.MosekSolver(LOG=0))
quad_mip_solvers = copy(quad_solvers)
# Solvers that take SOC in quadratic form
quad_soc_solvers = Any[]
grb && push!(quad_soc_solvers, Gurobi.GurobiSolver(QCPDual=1,OutputFlag=0))
cpx && push!(quad_soc_solvers, CPLEX.CplexSolver(CPX_PARAM_SCRIND=0))
xpr && push!(quad_soc_solvers, Xpress.XpressSolver(OUTPUTLOG=0, FEASTOL = 1e-9, BARPRIMALSTOP = 1e-9, BARGAPSTOP = 1e-9, BARDUALSTOP = 1e-9))
#osl && push!(quad_solvers, CoinOptServices.OsilSolver(CoinOptServices.OSOption("sb","yes",solver="ipopt")))
soc_solvers = copy(quad_solvers)
ipt && push!(quad_solvers, Ipopt.IpoptSolver(print_level=0))
Expand Down

0 comments on commit 3e66b34

Please sign in to comment.