Convex.jl for SDP

Specific Domains Optimization (Mathematical)
,
1

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

image
image1044×485 45.6 KB

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()
   @ Main ~/.julia/dev/DDControl/bl_ctrl_quad_stab.jl:29
 [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

[1] Khlebnikov, Mikhail V. “Quadratic stabilization of bilinear control systems.” Automation and Remote Control 77 (2016): 980-991.

2

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
3

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