I am using SOS programming to find a Lyapunov function, obeying some natural symmetry. The time derivative of V also obeys said symmetry. The feasibility of the Lyapunov function is dependent on the largeness of a parameter specified at the beginning of the problem which determines the dynamics of the differential equation. I will outline a sort of pseudocode here.
Truncated system:
@polyvar a[1:6]
Vbasis = (some vector of symmetry-invariant polynomial terms in a)
V_coeffs = @variable(model, [1:length(Vbasis)])
V = sum(V_coeffs.*Vbasis)
dVdt = f.*differentiate(V,a)
model = SOSModel(Mosek.Optimizer)
con_Vsos = @constraint(model, V - epsilon*dot(a,a) >= 0, symmetry=pattern)
minusdVdtsos = -dVdt - epsilon*dot(a,a)
con_dv = @constraint(model, minusdVdtsos >= 0, symmetry=pattern)
optimize!(model)
This works beautifully, the symmetry is imposed and outperforms the setting where we enforce sparsity=SignSymmetry() instead. Now, I add another term to the Lyapunov function V for which I only have estimates on dV/dt.
@polyvar a[1:6], q
Vbasis = (some vector of symmetry-invariant polynomial terms in a, q^2)
V_coeffs = @variable(model, [1:length(Vbasis)])
V = V_coeffs.*Vbasis
sbasis = (some vector of polynomial terms in a, q^2)
s_coeffs = @variable(model, [1:length(sbasis)])
s = sum(s_coeffs.*sbasis
bound = constant*s*q^2 - (positive term)
Mis = differentiate(V,a) - 2*dV/d(q^2)*a
model = SOSModel(Mosek.Optimizer)
# These constraints enforce the fact that |Mis[1]| <= s; we do this for the length of Mis (6) but I only put in one constraint for simplicity
@constraint(model, s + Mis[1]>=0)
@constraint(model, s - Mis[1] >=0)
dVdt = f.*differentiate(V,a) + bound
con_Vsos = @constraint(model, V - epsilon*dot([a;q],[a;q]) >= 0, symmetry=pattern)
minusdVdtsos = -dVdt - epsilon*dot([a;q],[a;q])
con_dv = @constraint(model, minusdVdtsos >= 0, symmetry=pattern)
optimize!(model)
If I do not enforce the natural symmetry on V, the program runs beautifully, and indeed the threshold parameter I am using to determine feasibility is equivalent between this Julia program and the old Matlab program I am running. Moreover, if I enforce instead sparsity=SignSymmetry() on any one of the constraints (turning them on and off like circuit breakers), it also preserves the feasibility result below the threshold parameter. However, the minute I enforce the natural symmetry, Mosek returns an infeasibility result unless my parameter is way below the threshold parameter. There are similar symmetries that can also be imposed on the bounds, but those also behave badly.
Here are multiple checks that have been made so far:
- The natural symmetry with q is the same as the natural symmetry without q, we do not enforce any symmetries on q so the pattern generation is identical except for the q term. I have verified that
- the homomorphism is indeed a homomorphism
- both V and minusdVdtsos are invariant under the defined homomorphism - I have ensured that any coefficients on any variables for the program are all above a certain threshold tolerance so there is no extra noise being sent to the Mosek solver. (There was actually a non-trivial observation here that caused this to be not as straightforward as I expected, and I may post about it in a different thread to ask for an error request for other users at some point, hopefully it’s a useful observation.)
- The size of the bases for the constraints in both the unsymmetrized (monomial basis) and symmetrized (nonmonomial basis) problems are of the same length, so the problem has many fewer constraints and variables that the original non-symmetrized problem. (Easy to be seen in the Mosek output as well.)
- I have minimized the number of free variables for the SOS program; for example, if V(a, q) = c_1(a1^2 + a2^2) + c2q^2, the free variables are also written this way.
My best intuition is there is some way in which JuMP is formulating the problem to hand to Mosek that causes Mosek to behave much worse numerically than it should. Supporting this intuition, I will sometimes get warnings from Mosek that the matrix ‘A’ is sparse. Indeed, this error also sometimes ocurs when I enforce the restriction sparsity=SignSymmetry(), but it does not strongly hamper the feasibility results like my handwritten symmetry. I was concerned about the basis size for the original problem vs the symmetrized problem, if one was longer than the other, but they were the same size so I was not worried about the fact that the block symmetry does not enforce Newton polytope reduction. Any time you have to peek at this would be greatly appreciated. Thanks.
