Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

don't automatically transform quadratic constraints into conic form #999

Merged
merged 1 commit into from
Mar 28, 2017
Merged

don't automatically transform quadratic constraints into conic form #999

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
123 changes: 4 additions & 119 deletions src/solvers.jl
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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