jump-dev / SumOfSquares.jl

Sum of Squares Programming for Julia
Other
118 stars 24 forks source link

How can we add some constraints? #19

Closed jonaprieto closed 5 years ago

jonaprieto commented 7 years ago

Hi, let's consider the following problem:

x1' = a - b*(x1 + 1)*x2 - a*x2
x2' = b * (x1+1) *x2 - c*x2

We want to find a Lyapunov function as you did in sosdemo2.jl file, so we only modified the vector field, and also the first inequality. Then, after these modifications, the file looks like

...
@polyvar x[1:2]
a = 2
b = 1
c = 1
f = [  a - b*(x[1] + 1)*x[2] - a*x[2],
        b * (x[1]+1) *x[2] - c*x[2] ]
m = Model(solver = solver)
Z = x.^2
@polyvariable m V Z
@polyconstraint m V >= 1e-16 * sum(x.^2)
P = dot(differentiate(V, x, f)
@polyconstraint m P <= 0
status = solve(m)
...

where f is the problem described above (vector field) and the solver iscsdp or scs.

As you can see, I defined the values for a, b, and c variables. Nevertheless, I'd like to add constraints over these variables but we didn't figure out how:

1 <= a <= 5
0 <= b <= 2
1  <= c <= 3

and it will be great also add contraints to the domain

 -1 <= x[1] <= 0
0   <= x[2] <= 1

Can you help us, please?

CC'ing @polislizarralde

blegat commented 7 years ago

Hi, thanks for reaching out to us. Could you be more precise about what you are trying to achieve ? You want to find a common Lyapunov function valid for all the systems for which 1 <= a <= 5, 0 <= b <= 2 and 1 <= c <= 3 ? That is an infinite number of systems so it seems tricky. However, since f depends affinely on a, b and c:

f = a * [1 - x[2], 0] + b * [(x[1] + 1)*x[2], (x[1]+1) *x[2]] + c * [0, -x[2]] = a f1 + b f2 + c f3

then you need

p(x) = a⟨∇V, f1⟩ + b⟨∇V, f2⟩ + c⟨∇V, f3⟩

to be positive for all (a,b,c) ∈ [1,5] x [0,2] x [1,3]. Note that the set of possible values for (a,b,c) is a polytope. Any point of a polytope is convex combination of the extreme points so you know that (a,b,c) is a convex combination of (1,0,1), (1,0,3), (1,2,1), (1,2, 3), (5,0,1), (5,0,3), (5,2,1), (5,2,3). Therefore p(x) is also the convex combination of the corresponding 8 polynomials. Long story short, if the Lyapunov works for the 8 dynamics then it works for all the possible values of (a,b,c). I have created this package for computing Lyapunov for switched systems.

For your domain constraints you can use the domain keyword like so

@polyconstraint(..., domain = -1 <= x[1] && x[1] <= 0 && 0 <= x[2] && x[2] <= 1)

See these examples for the domain keyword.

blegat commented 7 years ago

If a, b and c are constant/cannot vary in time, you do not need the same Lyapunov for every a, b and c. In that case you can declare them as polyvar and use them in V(x). For example:

@polyvar a b c
@polyvariable m V monomials([x; a; b; c], 0:5)

a polynomial V with all the possible monomials containing x[1], x[2], a, b or c with a degree between 0 and 5. Then of course you need to add the constraints on a, b and in the domain when building the constraints.

blegat commented 7 years ago

@jonaprieto Any update ?

blegat commented 5 years ago

Closing for inactivity, please reopen if you have any update