Transcript No Slide Title

Stability Region Analysis using composite Lyapunov
functions and bilinear SOS programming
Support from AFOSR
FA9550-05-1-0266, April 05-November 06
Authors
Weehong Tan, Andy Packard
Mechanical Engineering, UC Berkeley
Acknowledgements
Thanks to Ufuk Topcu, Gary Balas and Pete Seiler; PENOPT
Website
http://jagger.me.berkeley.edu/~pack/certify
Copyright 2006, Packard, Tan, Wheeler, Seiler and Balas. This work is licensed under the Creative Commons Attribution-ShareAlike License.
To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/2.0/ or send a letter to Creative Commons, 559 Nathan Abbott
Way, Stanford, California 94305, USA.
Quantitative Nonlinear Analysis
Initial focus
– Region of attraction estimation
– Attractive invariant sets
– L2  L2 induced norms
– L2  L induced norms
for
– finite-dimensional nonlinear systems, with
• polynomial vector fields
• parameter uncertainty (also polynomial)
Main Tools:
– Lyapunov/HJI formulation
– Sum-of-squares proofs to ensure nonnegativity and set containment
– Semidefinite programming (SDP), Bilinear Matrix Inequalities
• Optimization interface: YALMIP and SOSTOOLS
• SDP solvers: Sedumi
• BMIs: using PENBMI (academic license from www.penopt.com)
Estimating Region of Attraction
Dynamics, equilibrium point
x  f ( x),
f (x)  0
V 1
User-defined function p, whose
sub-level sets are to be in
region-of-attraction
P : x : p( x)     ROA x
p3
p2
p 1
x
By choice of positive-definite V,
maximize  so that
x : p( x)     x:V(x) 1
 dV

x : x  x ,V ( x)  1   x: f  0
 dx

dV
f 0
dx
Convexity of Analysis
In a global stability analysis, the certifying Lyapunov
functions
V : V  0
with V  f  0
are themselves a convex set.
In local analysis, the condition holds on sublevel sets
V : V  0
with V  f  0 on V  1
This set of certifying Lyapunov functions is not convex.
Example: f ( x)   x
V1 ( x)  16 x 2  19.95 x 3  6.4 x 4
V2 ( x)  0.1x 2
Vc ( x)  0.58  V1 ( x)  0.42  V2 ( x)
Estimating Region of Attraction
Dynamics, equilibrium point
x  f ( x),
f (x)  0
V 1
User-defined function p, whose
sub-level sets are to be in
region-of-attraction
P : x : p( x)     ROA x
p3
p2
p 1
x
By choice of positive-definite V,
maximize  so that
x : p( x)     x:V(x) 1
 dV

x : x  x ,V ( x)  1   x: f  0
 dx

dV
f 0
dx
Sum-of-Squares
Sum-of-squares decompositions will be the main tool to
decide set containment conditions, and certify nonnegativity.
A polynomial f, in n real-variables is a sum-of-squares if it
can be expressed as a sum-of-squares of other polys,
p
f   g 2j
j 1
Notation
 n set of all sum-of-square polynomials in n variables
Pn set of all polynomials in n variables
Sum-of-Squares Decomposition
For a polynomial f, in n real-variables, and of degree 2d
f  n

M 0 such that f  z T Mz
where z  [1, x1 , x2 ,, xn , x1 x2 ,, xnd ]T .
The entries of z are not algebraically independent
– e.g. x12x22 = (x1x2)2
– M is not unique (for a specified f)
The set of matrices, M, which yield f, is an affine subspace
– one particular + all homogeneous
– Particular solution depends on f
– all homogeneous solutions depend only on n & d.
Searching this affine subspace for a p.s.d element is a SDP…
Sum-of-Squares as SDP
For a polynomial f, in n real-variables, and of degree 2d
f  n

q
  R q such that M 0   i M i  0
i 1
Semidefinite program: feasibility
Each Mi is s×s, where
n  d 
s

 d 
2

n

d
1 
  n  d    n  2d 
q  
 
  

2  d   d   2d 
Using the Newton polytope method, both s and q can
often be reduced, depending on the terms present in f.
(s,q) dependence on n and 2d
n  d 

s
 d 


2d
n
 n  d  2  n  d   n  2d 
1 
 
  

q
2  d   d   2d 


2
4
6
8
2
3
0
6
6
10
27
15
75
3
4
0
10
20
20
126
35
465
4
5
0
15
50
35
420
70
1990
6
7
0
28
196
84
2646
210
19152
8
9
0
45
540 165 10692
495
109890
10
11
0
66 1210 286 33033 1001
457743
Synthesizing Sum-of-Squares as SDP
Given: polynomials f 0 , f1 ,, f m
Decide if an affine combination of them can be made a
sum-of-squares.
This is also an SDP.
m
  R m with f 0   k f k   n
k 1

qm
  R q  m with M 0   i M i  0
i 1
Synthesizing Sum-of-Squares as Bilinear SDP
Given: polynomials
f 0 , f1 , , f m
g 0 , g1 ,  , g m
h0 , h1 ,  , hm
A problem that will arise in this talk is: find
such that
  R m ,  R m
m
m



f 0   k f k   g 0   k g k  h0   k hk    n
k 1
k 1
k 1



m
This is a nonconvex SDP, namely a bilinear matrix inequality
Psatz
Given: polynomials f1 , , f m , g1 , , g p , h1 , , hq
 x : f1 ( x)  0, , f m ( x)  0, 
Goal: Decide if


 :  g1 ( x)  0, , g p ( x)  0,
the set is empty.
 h1 ( x)  0, , hq ( x)  0 


Φ is empty if and only if
f  P  f1 ,, f m , g  M g1 ,, g p , h  I h1 ,, hq 
such that f  g 2  h  0
 g :   Z 
P  f ,, f    s b : s  , b  M  f ,, f 
I h ,, h    h p : p R 
M g1 ,, g p  
1
1
i
i
m
q
p
i i
k
k
i
i
k
1
m
Region of Attraction
V 1
p3
p2
p 1
By choice of positive-definite V,
maximize  so that
x : p( x)     x:V(x) 1
x : x  x,V ( x)  1  x:V  f
 0
dV
f 0
dx
max  over V , V ( x )  0, s6 , s8 , s9   n
“small” positive
subject to
definite functions
V  l1   n
BMIs
   p s6  V  1   n
 1  V s8  s9V  f  l2    n
Simple Psatz:
PENBMI from
PENOPT
x
Products of
decision variables
Sanity check
For a positive definite matrix B,
x   x  xBx xROA  x : xBx  1
Proof: V x  : xBx
0
 V  V V  1.
nth order system
cubic vector field
known ROA
Consider p.d. quadratic shape factor p( x)  xT Rx
The best obtainable result is the “largest” value max such that
x : xRx     x:xBx  1
That containment easy to characterize:
 
 12
max  max R BR
 12

1
Questions:
Yes
– Can the formulation we wrote yield this?
– Can the BMI solver find this solution? Basically, Yes
100’s of random examples, n=2-8; 3 restarts
of PENBMI, always successful
Example: Van der Pol: ROA
Classical 2-d system
x1   x2
x2  x1  x12  1 x2


Except for nV=4 case, the
results are comparable to
Papachristodoulou (2005) and
Wloszek (2003)
V # decvars
2
13
4
57
6
166
Features:
– Unstable limit cycle around origin
– One equilibrium point: stable, at origin
– Here, we use an elliptical shape factor
Region of Attraction: pointwise-max
If V1 and V2 are positive definite, and
V1 ( x)  1
dV1
f ( x)  0 on
V2 ( x)  V1 ( x)
dx
and dV2
V2 ( x)  1
f ( x)  0 on
V1 ( x)  V2 ( x)
dx
V2  1
V1  1
V 1
Then
V ( x) : max V1 ( x), V2 ( x)
proves asymptotic stability of
x  f (x)
on
x : V ( x)  1
Region of Attraction with pointwise-max
Use Psatz to get a sufficient condition for
x : p( x)     x:V(x) 1
x : x  0,V ( x)  1  x:V  f
using V of the form
 0
V ( x) : max V1 ( x), V2 ( x), ,Vq ( x)
i
max  over Vi , Vi (0)  0, s6i , s8i , s9i , s0ij   n
subject to
Vi  l1   n
   p s6i  Vi  1   n
 1  Vi s8i  s9i Vi  f  l2  
 s V  V  
q
j 1, j  i
0 ij
i
j
n
ROA with Pointwise-Max Lyapunov functions
Vi q  1 q  2
2 13
38
4
57 120
6 166 338
Original
(single
V)i)
Composite
(2 V
Pointwise max of 6th degree V1,V2
I1, A
V1  1, V2  1  V1
V1  0
V1  1
V2  1

V1  1, V2  1  V1
V  0
I1, B
1
Is V2 , by itself, a decent Lyapunov function?
– Sub-Level set looks similar to result,
  x : V2 ( x)  1
– But, derivative on sublevel set is not negative
V2  0
V2  0
Different shape factor
Nearly the
same results.
ROA: 3rd order example
Example (from Davison, Kurak):
x1   x2
x2   x3
 12.5  8.1 3.0
p( x)  xT  8.1 20.8  8.5 x


 3.0  8.5 13.4
x3  0.915 x1  1  0.915 x12 x2  x3

Solutions diverge from
these initial
conditions,i.e these
initial conditions are
not in the ROA

Problems, difficulties, risks
Dimensionality:
– For general problems, it seems unlikely to move beyond cubic
vector fields and (pointwise-max) quadratic V. These result in
“tolerable” SDPs for state dimension < 15.
– Theory may lead to reduced complexity in specific instances of
problems (sparsity, Newton polytope reduction, symmetries)
Solvers (SDP): numerical accuracy, conditioning
Connecting the Lyapunov-type questions to MilSpec-type
measures
– Decay rates
– Damping ratios
– Oscillation frequencies
– Time-to-double
BMI nature of local analysis