Did you have a look at JuMP, especially the Section Semidefinite Constraints or the Example Example minimal Ellipse?
I have no experience so far in LMI’s but I think this would be a good starting point.
For example solving the Lyapunov Inequality for a continuous system could look like this:
using JuMP
using MosekTools
using LinearAlgebra
A = [-0.5 1; 0 -0.3]; n = size(A,1); Q = [1 0; 0 1]
model = Model(with_optimizer(Mosek.Optimizer))
@variable(model, P[i=1:n, j=1:n], PSD)
@SDconstraint(model, P ⪰ 0)
@SDconstraint(model, A'*P + P*A ⪯ -Q)
@SDconstraint(model, A'*P + P*A ⪰ -Q)
optimize!(model)
P = value.(P)
# Test result:
A'*P + P*A ≈ -Q