# Convex.jl gives solution but constraints are not satisfied

**URL:** <https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428>\
**Category:** Optimization (Mathematical)\
**Tags:** optimization, convex-optimization, convex\
**Created:** [June 16, 2023, 12:03am UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428 "2023-06-16T00:03:56Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [June 16, 2023, 12:03am UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/1 "2023-06-16T00:03:56Z")

</div>

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:

```julia
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];
    A = [Ad,Ad];
    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
    add_constraint!(Q1,Q1 in :SDP); add_constraint!(Q2,Q2 in :SDP);

    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

```

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [June 17, 2023, 9:19am UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/2 "2023-06-17T09:19:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [June 17, 2023, 10:03am UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/3 "2023-06-17T10:03:47Z")

</div>

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:

```julia
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.

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 19, 2023, 12:12am UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/4 "2023-06-19T00:12:47Z")

</div>

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

```plaintext
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:

```plaintext
using JuMP, LinearAlgebra, SCS

function feedback_jump(ρ)
    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]
    A = [Ad, Ad]
    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)

```

---

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [June 20, 2023, 8:16pm UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/5 "2023-06-20T20:16:59Z")

</div>

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 20, 2023, 8:53pm UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/6 "2023-06-20T20:53:35Z")

</div>

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

```julia
E = eigen(Matrix(si3))
E.values .= max.(0.0, E.values)
si3_new = Matrix(E)

```

---

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [June 20, 2023, 9:55pm UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/7 "2023-06-20T21:55:23Z")

</div>

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`

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  
]

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [June 20, 2023, 10:29pm UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/8 "2023-06-20T22:29:57Z")

</div>

> [@shaurya](#):
>
> As in, it showing the problem is solved when one of the constraints is not satisfied when checked manually later.

SCS is a first-order solver. It uses these termination conditions: [Algorithm — SCS 3.2.4 documentation](https://www.cvxgrp.org/scs/algorithm/index.html#optimality). 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](https://www.cvxgrp.org/scs/api/settings.html).

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](https://jump.dev/JuMP.jl/stable/packages/Clarabel/).

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

---

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [June 23, 2023, 8:08pm UTC](https://discourse.julialang.org/t/convex-jl-gives-solution-but-constraints-are-not-satisfied/100428/9 "2023-06-23T20:08:50Z")

</div>

Thanks! That helps! 🙂