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];
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
