# Convex.jl gives solution but constraints are not satisfied

Hi,

I’ve been using convex.jl for some of my problems but I feel stuck in one SDP implementation. When I run the following code, the solver says that it’s feasible and reaches an optimal point, however, when I check the last constraint ( commented as equation 4 in the code) if it is Positive semi definite, it returns false. This is surprising and I’m not sure how to proceed about it. Thanks for your help, I’ve attached the code below:

``````using Convex, LinearAlgebra, SCS

function feedback(ρ)

# data matrices
Ad  = [ 1.0         0.0   1.0         0.00114404
-1.49734e-9  1.0  -0.00114404  0.999999
3.92645e-6  0.0   0.999999    0.00228807
-4.492e-9    0.0  -0.00228807  0.999997];
Bd  = [ 0.5          0.000381345
-0.000381345  0.5
1.0          0.00114404
-0.00114404   0.999999];
P   = [ 0.00383029  0.995214
0.995214    0.00478556];
B       = [Bd,zeros(4,2)];
ic      = [-1.5;-0.5;5e-3;0];

# defining optimization model
Q1 = Variable(4,4); Q2= Variable(4,4);
Q= [Q1,Q2];  δ = Variable();
Y1 = Variable(2,4); Y2= Variable(2,4);
Y= [Y1,Y2];
prob = minimize(δ);

C = [Matrix(0.01I,4,4);zeros(2,4)];
D = [zeros(4,2);Matrix(1.0I,2,2)];

# setting PSD constraints

for i in 1:2
AQBY    = A[i]*Q[i] + B[i]*Y[i];
Γ       = [sqrt(P[i,j])*Matrix(1.0I,4,4)     for j in 1:2];
Γ   =  reduce(hcat,Γ);
# equation 1
FI      = [1        ic'
ic      Q[i] ];
prob.constraints += (FI in :SDP);

# equation 2
Qₑ = [Q1 zeros(4,4); zeros(4,4)     Q2];
SI = [          Q[i]        AQBY'*Γ                                     Q[i]*C'                                             Y[i]'*D'
Γ'*AQBY         Qₑ                                          zeros(size(Qₑ,1),size(Q[i]*C',2))                   zeros(size(Qₑ,1),size(Y[i]'*D',2))
C*Q[i]      zeros(size(C*Q[i],1),size(Qₑ,2))        δ*Matrix(1.0I,size(C*Q[i][:,:],1),size(Q[i]*C',2))       zeros(size(C*Q[i],1),size(Y[i]'*D',2))
D*Y[i]      zeros(size(D*Y[i],1),size(Qₑ,2))        zeros(size(D*Y[i],1),size(Q[i]*C',2))           δ*Matrix(1.0I,size(D*Y[i],1),size(Y[i]'*D',2))  ];
prob.constraints += (SI in :SDP);

# equation 3
for j in 1:2
SI2 = [     Q[i][:,:]    AQBY'
AQBY        Q[j][:,:]];
prob.constraints += (SI2 in :SDP);
end

end

# equation 4
G=Matrix(1.0I,2,2);
SI3 = [ ρ^2*Matrix(1.0I,2,2)       G*Y1
Y1'*G'                      Q1'   ];
prob.constraints += (SI3 in :SDP);

Convex.solve!(prob, SCS.Optimizer; silent_solver = false)

# checking the constraint
Y_sol = evaluate(Y1)
Q_sol = evaluate(Q1);
si3 = [ ρ^2*Matrix(1.0I,2,2)       G*Y_sol
Y_sol'*G'                 Q_sol'   ];

return prob,isposdef(si3),si3

end

``````

It’s a wild guess, but maybe the matrix is at the limit between SDP and not SDP, and numerical errors (or tolerances) cause the SDP test to fail.

1 Like

One way to check cvanaret’s guess would be to use `eigen` to get all eigenvalues and see if they are positive and real enough.

For example:

``````using LinearAlgebra

si3 = rand(3,3)

# check if eigenvalues are real enough:
maximum(abs.(imag.(eigen(si3).values)))
# above should be a small number. works surely for si3*si3':
maximum(abs.(imag.(eigen(si3*si3').values)))

# check if eigenvalues are non-negative
# (i.e. minimum should be positive):
minimum(real.(eigen(si3*si3').values))
# should be okay for si3*si3' as it is positive semidefinite.
``````
2 Likes

The issue is that the `si3` matrix has positive eigen values but is not symmetric.

``````julia> prob, _, si3 = feedback(1.0);
------------------------------------------------------------------
SCS v3.2.3 - Splitting Conic Solver
(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 50, constraints m: 1527
cones: 	  z: primal zero / dual free vars: 712
s: psd vars: 815, ssize: 11
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
max_iters: 100000, normalize: 1, rho_x: 1.00e-06
acceleration_lookback: 10, acceleration_interval: 10
lin-sys:  sparse-direct-amd-qdldl
nnz(A): 1092, nnz(P): 0
------------------------------------------------------------------
iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
0| 1.81e+00  1.00e+00  9.38e-01  1.76e-01  1.00e-01  8.13e-04
250| 1.66e-03  2.68e-04  1.55e-03  3.36e-03  1.00e-01  6.00e-02
500| 8.77e-04  1.81e-04  1.77e-03  3.48e-03  1.00e-01  1.17e-01
750| 8.05e-04  1.98e-04  1.97e-03  3.39e-03  1.00e-01  1.73e-01
1000| 7.51e-04  2.15e-04  2.12e-03  3.32e-03  1.00e-01  2.37e-01
1250| 7.15e-04  2.28e-04  2.19e-03  3.26e-03  1.00e-01  3.04e-01
1500| 6.89e-04  2.36e-04  2.20e-03  3.21e-03  1.00e-01  3.63e-01
1750| 6.69e-04  2.39e-04  2.17e-03  3.19e-03  1.00e-01  4.22e-01
2000| 6.53e-04  2.39e-04  2.09e-03  3.19e-03  1.00e-01  4.84e-01
2250| 6.40e-04  2.34e-04  1.97e-03  3.20e-03  1.00e-01  5.45e-01
2500| 6.28e-04  2.24e-04  1.82e-03  3.24e-03  1.00e-01  6.03e-01
2750| 8.32e-03  1.42e-03  1.85e-04  4.68e-03  1.00e-01  6.67e-01
3000| 6.09e-04  1.94e-04  1.43e-03  3.38e-03  1.00e-01  7.24e-01
3250| 6.00e-04  1.74e-04  1.20e-03  3.48e-03  1.00e-01  7.82e-01
3500| 5.92e-04  1.52e-04  9.56e-04  3.60e-03  1.00e-01  8.40e-01
3750| 5.84e-04  1.28e-04  7.02e-04  3.72e-03  1.00e-01  9.07e-01
4000| 5.77e-04  1.02e-04  4.44e-04  3.85e-03  1.00e-01  9.63e-01
4250| 5.71e-04  7.74e-05  1.92e-04  3.99e-03  1.00e-01  1.02e+00
4350| 5.68e-04  6.77e-05  9.45e-05  4.04e-03  1.00e-01  1.04e+00
------------------------------------------------------------------
status:  solved
timings: total: 1.04e+00s = setup: 7.03e-04s + solve: 1.04e+00s
lin-sys: 8.73e-02s, cones: 9.29e-01s, accel: 5.06e-03s
------------------------------------------------------------------
objective = 0.004039
------------------------------------------------------------------

julia> eigen(si3)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
6-element Vector{Float64}:
0.003652162698009067
0.017117415129939587
0.031492005016552194
1.0000144031866585
1.000018360687806
4.21695046889853
vectors:
6×6 Matrix{Float64}:
-0.00100689  -0.0040127   -0.00125398    0.0902535     0.995909     -9.83063e-5
0.00358184  -0.00117339   0.000556168  -0.995912      0.0902533    -1.05646e-5
-0.0321066    0.0144825    0.327017      9.05741e-5    0.000336223  -0.944362
0.0930193   -0.204946    -0.918735     -0.000106781  -0.00191088   -0.324449
-0.319671    -0.932906     0.157682     -0.00031196   -0.00385018    0.0511625
0.942397    -0.295729     0.155309      0.00378947   -0.000384923   0.017206

julia> si3 .- si3'
6×6 Matrix{Float64}:
0.0  0.0  0.0           0.0           0.0          0.0
0.0  0.0  0.0           0.0           0.0          0.0
0.0  0.0  0.0          -2.06907e-10  -6.27553e-9  -2.20299e-9
0.0  0.0  2.06907e-10   0.0          -1.12485e-9  -3.91686e-10
0.0  0.0  6.27553e-9    1.12485e-9    0.0          2.93937e-9
0.0  0.0  2.20299e-9    3.91686e-10  -2.93937e-9   0.0
``````

If you use JuMP, you can explicitly declare matrices as symmetric:

``````using JuMP, LinearAlgebra, SCS

function feedback_jump(ρ)
1.0 0.0 1.0 0.00114404
-1.49734e-9 1.0 -0.00114404 0.999999
3.92645e-6 0.0 0.999999 0.00228807
-4.492e-9 0.0 -0.00228807 0.999997
]
Bd = [0.5 0.000381345; -0.000381345 0.5; 1.0 0.00114404; -0.00114404 0.999999]
P = [0.00383029 0.995214; 0.995214 0.00478556]
B = [Bd, zeros(4, 2)]
ic = [-1.5; -0.5; 5e-3; 0]
model = Model(SCS.Optimizer)
@variable(model, Q1[1:4, 1:4], PSD)
@variable(model, Q2[1:4, 1:4], PSD)
@variable(model, Y1[1:2, 1:4])
@variable(model, Y2[1:2, 1:4])
@variable(model, δ)
@objective(model, Min, δ)
Q, Y = [Q1, Q2], [Y1, Y2]
C, D = [0.01I(4); zeros(2, 4)], [zeros(4, 2); I(2)]
for i in 1:2
AQBY = A[i] * Q[i] + B[i] * Y[i];
Γ = reduce(hcat,  [sqrt(P[i, j]) * I(4) for j in 1:2])
Qₑ = [Q1 zeros(4, 4); zeros(4, 4) Q2]
SI = [
Q[i]    AQBY'*Γ     Q[i]*C'     Y[i]'*D'
Γ'*AQBY Qₑ          zeros(8, 6) zeros(8, 6)
C*Q[i]  zeros(6, 8) δ*I(6)      zeros(6, 6)
D*Y[i]  zeros(6, 8) zeros(6, 6) δ*I(6)
]
@constraints(model, begin
Symmetric([1 ic'; ic Q[i]]) >= 0, PSDCone()
Symmetric(SI) >= 0, PSDCone()
[j=1:2], Symmetric([Q[i] AQBY'; AQBY Q[j]]) >= 0, PSDCone()
end)
end
G = I(2)
@constraint(model, Symmetric([ρ^2*G G*Y1; Y1'*G Q1']) >= 0, PSDCone())
optimize!(model)
Y_sol = value.(Y1)
Q_sol = value.(Q1)
si3 = [ρ^2*G (G*Y_sol); Y_sol'*G' Q_sol']
return model, isposdef(si3), si3
end

feedback_jump(1.0)
``````
1 Like

Thanks for your reply! I tried this and it works when I run for `feedback(1.0)`

However, it doesn’t give the desired results for `feedback(1e-3)`

I did this after running your JuMP function:

``````model,a,si3 = feedback_jump(1e-3);
solution_summary(model)
``````

The variable a shows that the matrix si3 has negative eigenvalues, although in the solution summary it shows the problem is solved with a feasible point! I don’t quite understand how that’s happening given it’s a constraint in the optimization.

Solvers don’t use exact arithmetic, so they satisfy constraints only up to a certain tolerance. I assume you’re getting small negative eigen values like `1e-6`?

You have a couple of options:

• Use the result as-is
• Re-scale the data. The problem is likely the `1e-9` magnitude coefficients in `Ad`. Solvers generally like all coefficients to be between `1e-4` and `1e4`.
• Project out the negative eigen values
``````E = eigen(Matrix(si3))
E.values .= max.(0.0, E.values)
si3_new = Matrix(E)
``````
1 Like

I set the 1e-9 magnitude coefficients to zero. That changes the problem slightly but it no longer has any coefficients beyond 1e-4 so should give better results.

If I run the same code again, I still have negative eigenvalues, with one of the eigenvalues = -0.00046. I’m not sure if this order suffices, but I’m still confused why this is happening for a convex problem. As in, it showing the problem is solved when one of the constraints is not satisfied when checked manually later.

For reference, I used the following `Ad`

1.0 0.0 1.0 0.00114404
0 1.0 -0.00114404 0.999999
0 0.0 0.999999 0.00228807
0 0.0 -0.00228807 0.999997
]

SCS is a first-order solver. It uses these termination conditions: Algorithm — SCS 3.2.4 documentation. It does not check that the constraints hold exactly, but that they hold to within some tolerance.

The default tolerances are quite loose (`eps_abs` and `eps_rel` are both `1e-4`): Settings — SCS 3.2.4 documentation.

You could try setting tighter tolerances with `set_attribute(model, "eps_abs", 1e-6)`, but this often causes SCS to not converge.

You could also try a different solver, like Clarabel: oxfordcontrol/Clarabel.jl · JuMP.

Just swap `using SCS` for `using Clarabel`, and `SCS.Optimizer` for `Clarabel.Optimizer`.

2 Likes

Thanks! That helps!

1 Like