Disclaimer: This is a repost/extension of this github issue.
The “output variables” of my problem are the entries in some specific place of a SDP matrix (P[1,4:5]) in the code below).
When I directly optimize an objective function over these indices of the matrix, I get the correct result.
However, when I then introduce new variables y[1:2] .== P[1,4:5] and re-solve the same objective, Mosek returns a different result.
To reproduce the error:
using JuMP, MosekTools, SCS
# setup the model and constraints
model = Model(Mosek.Optimizer)
@variable(model, P[1:5,1:5], PSD)
@constraint(model, diag(P[2:3, 2:3]) .- (-1. .+ 1) .* P[1,2:3] .+ (-1.) .* 1. .<= 0)
@constraint(model, P[1,1] .== 1.) W = [1. 1; 1 -1]
@constraint(model, P[1,4:5] .>= 0.0)
@constraint(model, P[1,4:5] .>= W*P[1,2:3])
@constraint(model, diag(P[4:5,4:5] .- W*P[2:3,4:5]) .== 0);
# solve the original objective
c = [1.14, -4.5]
@objective(model, Max, c'*P[1,4:5])
optimize!(model)
objective_value(model) # returns ~2.28
# reformulate objective
@variable(model, y[1:2])
@constraint(model, y .== P[1,4:5])
@objective(model, Max, c'*y)
optimize!(model)
objective_value(model) # returns ~4.55 !!!
When using model = Model(SCS.Optimizer) instead, the objective value is ~2.28 both times.
I’m using Julia 1.11.1, JuMP v1.23.6, MosekTools v0.15.4 and SCS v2.0.2.