Fix creation of polynomials with grammonomials keyword

Author: blegat

src/constraint.jl

@@ -30,7 +30,6 @@ addpolynonnegativeconstraint{T<:VectorOfPolyType{false}}(m::JuMP.Model, P::Matri
 addpolynonnegativeconstraint{T<:VectorOfPolyType{true}}(m::JuMP.Model, P::Matrix{T}, domain::BasicSemialgebraicSet) = matconstraux(PolyVar{true}, m, P, domain)
 
 function addpolynonnegativeconstraint(m::JuMP.Model, p, domain::BasicSemialgebraicSet)
-    # TODO We might want to add this as a function in MultivariatePolynomials.jl
     mindeg, maxdeg = extdeg(p)
     for q in domain.p
         mindegq, maxdegq = extdeg(q)
@@ -44,6 +43,7 @@ function addpolynonnegativeconstraint(m::JuMP.Model, p, domain::BasicSemialgebra
         s = createnonnegativepoly(m, :Gram,  MonomialVector(vars(p), mind:maxd))
         p -= s*q
     end
+    # FIXME If p is a MatPolynomial, p.x will not be correct
     Z = getmonomialsforcertificate(p.x)
     slack = createnonnegativepoly(m, :Gram, Z)
     q = p - slack

src/variable.jl

@@ -6,7 +6,7 @@ function createpoly{C}(m::JuMP.Model, monotype::Symbol, x::MonomialVector{C})
     end
     gram = monotype == :Gram
     if gram
-        Z = dot(x, x).x
+        Z = (sum(x)^2).x
     else
         Z = x
     end
@@ -15,21 +15,27 @@ end
 createpoly(m::JuMP.Model, monotype::Symbol, x::Vector) = createpoly(m, monotype, MonomialVector(x))
 
 function createnonnegativepoly{C}(m::JuMP.Model, monotype::Symbol, x::MonomialVector{C})
-    if monotype == :Default
-        monotype = :Gram
-    end
-    gram = monotype == :Gram
-    if gram
-        Z = x
+    if isempty(x)
+        # Need MultivariatePolynomials v0.0.2
+        #zero(MatPolynomial{C, JuMP.Variable})
+        MatPolynomial(JuMP.Variable[], x)
     else
-        Z = getmonomialsforcertificate(x)
-    end
-    p = MatPolynomial{C, JuMP.Variable}((i, j) -> Variable(m, -Inf, Inf, :Cont), Z)
-    push!(m.varCones, (:SDP, p.Q[1].col:p.Q[end].col))
-    if !gram
-        # The coefficients of a monomial not in Z do not all have to be zero, only their sum
-        addpolyeqzeroconstraint(m, removemonomials(Polynomial(p), x))
+        if monotype == :Default
+            monotype = :Gram
+        end
+        gram = monotype == :Gram
+        if gram
+            Z = x
+        else
+            Z = getmonomialsforcertificate(x)
+        end
+        p = MatPolynomial{C, JuMP.Variable}((i, j) -> Variable(m, -Inf, Inf, :Cont), Z)
+        push!(m.varCones, (:SDP, p.Q[1].col:p.Q[end].col))
+        if !gram
+            # The coefficients of a monomial not in Z do not all have to be zero, only their sum
+            addpolyeqzeroconstraint(m, removemonomials(Polynomial(p), x))
+        end
+        p
     end
-    p
 end
 createnonnegativepoly(m::JuMP.Model, monotype::Symbol, x::Vector) = createnonnegativepoly(m, monotype, MonomialVector(x))

test/constraint.jl

@@ -1,4 +1,4 @@
-@testset "Non-symmetric matrix SOS constraint with $solver" for solver in sdp_solvers
+@testset "Non-symmetric matrix SOS constraint" begin
     @polyvar x
     m = SOSModel()
     @test_throws ArgumentError addpolynonnegativeconstraint(m, [1 x; -x 0], BasicSemialgebraicSet())

test/lyapunov.jl

@@ -0,0 +1,2 @@
+include("switchedsystems.jl")
+#include("sosdemo2.jl")

test/runtests.jl

@@ -4,16 +4,18 @@ using SumOfSquares
 using PolyJuMP
 using Base.Test
 
+include("variable.jl")
+include("constraint.jl")
+
 include("solvers.jl")
 
 include("certificate.jl")
-include("constraint.jl")
 
 include("motzkin.jl")
 
 # SOSTools demos
 include("sospoly.jl")
-include("sosdemo2.jl")
+include("lyapunov.jl")
 include("sosdemo3.jl")
 include("sosdemo4.jl")
 include("sosdemo5.jl")

test/switchedsystems.jl

@@ -0,0 +1,45 @@
+# See https://github.com/blegat/SwitchedSystems.jl
+
+# Example 2.8 of
+# P. Parrilo and A. Jadbabaie
+# Approximation of the joint spectral radius using sum of squares
+# Linear Algebra and its Applications, Elsevier, 2008, 428, 2385-2402
+# Inspired from the construction of:
+# Ando, T. and Shih, M.-h.
+# Simultaneous Contractibility.
+# SIAM Journal on Matrix Analysis & Applications, 1998, 19, 487
+# The JSR is √2
+
+@testset "[PJ08] Example 2.8 with $solver" for solver in sdp_solvers
+    @polyvar x[1:2]
+    A1 = [1 0; 1 0]
+    A2 = [0 1; 0 -1]
+    expected_ub = [√2, 1]
+    function testlyap(d, γ, expected_status)
+        m = SOSModel(solver = solver)
+        X = monomials(x, d)
+        @polyvariable(m, p, grammonomials = X)
+        # p strictly positive
+        q = sum(x.^(2*d))
+        @polyconstraint m p >= q
+        c1 = @polyconstraint m p(A1 * x, x) <= γ^(2*d) * p
+        c2 = @polyconstraint m p(A2 * x, x) <= γ^(2*d) * p
+
+        @test expected_status == solve(m; suppress_warnings=true)
+
+        if expected_status == :Infeasible
+            μ1 = getdual(c1)
+            μ2 = getdual(c2)
+
+            # The dual constraint should work on any polynomial.
+            # Let's test it with q
+            lhs = dot(μ1, q(A1 * x, x)) + dot(μ2, q(A2 * x, x))
+            rhs = dot(μ1, q) + dot(μ2, q)
+            @test 1e-6 * max(abs(lhs), abs(rhs)) + lhs >= rhs
+        end
+    end
+    testlyap(1, √2 - 1e-4, :Infeasible)
+    testlyap(1, √2 + 1e-3, :Optimal)
+    testlyap(2, 1 - 1e-3, :Infeasible)
+    testlyap(2, 1 + 1e-2, :Optimal)
+end

test/variable.jl

@@ -0,0 +1,6 @@
+@testset "Creating polynomial with empty MonomialVector" begin
+    x = MonomialVector{true}()
+    m = SOSModel()
+    @inferred createnonnegativepoly(m, :Gram, x)
+    @test createnonnegativepoly(m, :Gram, x) == 0
+end