# Convex.jl for SDP

I am trying to follow up [1, Example 1]. This involves semidefinite programming described in Corollary 1:

The following is what I’ve tried, and basically I think I need some specialized syntax for SDP.

# Trial

## Code

using Convex
using LinearAlgebra
using SCS

"""
[1, Example 1]
# Refs
[1] Khlebnikov, Mikhail V.
"Quadratic stabilization of bilinear control systems."
Automation and Remote Control 77 (2016): 980-991.
"""
function main()
A = [0 1;
1 1]
b = [0 1]'
D = [1 1;
-1 1]
n = 2
ϵ = 0.1  # TODO: variable?
μ = 0.0  # TODO: see remark of Corollary 1
P = Semidefinite(n)
y = Variable(n)
tmp = [(A*P + P*A' + b*y' + y*b' + ϵ*D*P*D') y;
y' -ϵ*I]
prob = minimize(
-eigmin(P),
isposdef(-tmp),
)
solve!(prob, SCS.Optimizer; silent_solver=true)
k = inv(evaluate(P)) * y
println(evaluate(P))
println(evaluate(k))
end


## Error message

julia> main()
ERROR: MethodError: no method matching isless(::Matrix{Any}, ::Int64)

Closest candidates are:
isless(::AbstractFloat, ::Real)
@ Base operators.jl:179
isless(::ForwardDiff.Dual{Tx}, ::Integer) where Tx
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:144
isless(::ForwardDiff.Dual{Tx}, ::Real) where Tx
@ ForwardDiff ~/.julia/packages/ForwardDiff/PcZ48/src/dual.jl:144
...

Stacktrace:
[1] <(x::Matrix{Any}, y::Int64)
@ Base ./operators.jl:343
[2] <=(x::Matrix{Any}, y::Int64)
@ Base ./operators.jl:392
[3] main()
[4] top-level scope
@ REPL[2]:1
[5] top-level scope
@ ~/.julia/packages/Infiltrator/LtFao/src/Infiltrator.jl:726


I followed this example provided by Convex.jl but somehow it fails.
Can anyone help me?

# Refs

To include SDP constraint, I needed to the following change, which is basically:

1. constructing matrix via vcat and hcat
2. appropriate clarification of matrix dimension (not \epsilon * I but epsilon for this case)
    tmp = vcat(hcat((A*P + P*A' + b*y' + y*b' + ϵ*D*P*D' + μ*P), y), hcat(y', -ϵ))
prob = minimize(
-eigmin(P),
isposdef(-tmp),
)


Note: this “seems” to be solved but the values are different from the reported ones in the reference.

1 Like

It should work if you do:

tmp = [
(A*P + P*A' + b*y' + y*b' + ϵ*D*P*D') y;
y' -ϵ*Matrix(I(1))
]


Convex doesn’t support UniformScaling from LinearAlgebra.

2 Likes