ROA calculations on the Van der Pol dynamics

Contents

This demo intends to show the basic problem setup and how one can manipulate the options in the ROA code. For more details on all the options see NACodeExplanation.pdf.

Default call

Form the vector field

pvar x1 x2;
x = [x1;x2];
x1dot = -x2;
x2dot = x1+(x1^2-1)*x2;
f = [x1dot; x2dot];

Get the default values of options to run the ROA code.

[roaconstr,opt,sys] = GetRoaOpts(f, x);

Call the wrapper which in turn calls RoaEst.m

outputs = wrapper(sys,[],roaconstr,opt);

Extract the solution

[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);

Certified beta

beta
beta =

    1.5138

Certifying V

V
 
V = 
  0.44698*x1^2 - 0.16958*x1*x2 
  + 0.39439*x2^2
 

Corresponding gamma

gamma
gamma =

    0.7712

Changing the options (p, degree of V, and type of iterations)

Default option in the ROA code uses

- simulation data and linearization for generating the initial V
- the iterations with no bisection
- p = x'*x
- quadratic V
- ...

Choice of p, the degree of V, and the type of iterations affects the way the options are set in GetRoaOpts.m (for example how the basis vectors for the multipliers are generated). Therefore, if these are to be set to non-default values, this should be done before the call to GetRoaOpts. For example,

Change the degree of multipliers

Bis.r1deg = 4;

V is created as V(x) = zV(x)'*alpha, where zV is a vector of monomials in x and alpha is a vector of free variables. The form V is therefore assigned by generating zV vector and passing into GetRoaOpts For example, the following generates zV as the vector of monomials of degree 2,3, and 4

zV = monomials(x,2:4)
 
zV = 
  [      x1^2 ]
  [     x1*x2 ]
  [      x2^2 ]
  [      x1^3 ]
  [   x1^2*x2 ]
  [   x1*x2^2 ]
  [      x2^3 ]
  [      x1^4 ]
  [   x1^3*x2 ]
  [ x1^2*x2^2 ]
  [   x1*x2^3 ]
  [      x2^4 ]
 

Now, call GetRoaOpts to get the correponding opts, roaconstr, etc.

[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p);

if any of zV or p is not changed from the default, just pass [].

Call the wrapper which in turn calls RoaEst.m

outputs = wrapper(sys,[],roaconstr,opt);

Extract the solution

[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);

Certified beta

beta
beta =

    2.0305

Certifying V

V
 
V = 
  0.031821*x1^4 + 0.026437*x1^3
  *x2 + 0.14315*x1^2*x2^2 
  - 0.10029*x1*x2^3 + 0.044507
  *x2^4 - 0.00012905*x1^3 
  + 0.00034277*x1^2*x2 
  - 0.00025672*x1*x2^2 
  - 0.00034085*x2^3 + 0.19464
  *x1^2 - 0.1828*x1*x2 + 0.11424
  *x2^2
 

Corresponding gamma

gamma
gamma =

    0.8084

Changing all other options

All other options can be changed after the call to GetRoaOpts (changing them will not affect other options). For example, changing option opt.coordoptim.IterStopTol to a smaller number will likely to improve the estimate because it will cause the coordinate-wise to stop after a larger number of iterations

opt.coordoptim.IterStopTol

opt.coordoptim.IterStopTol = 1e-5;
betaWithLargeStoppingTolerance = beta;

outputs = wrapper(sys,[],roaconstr,opt);
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
betaWithSmallStoppingTolerance = beta;

betaWithLargeStoppingTolerance
betaWithSmallStoppingTolerance
ans =

    0.0300


betaWithLargeStoppingTolerance =

    2.0305


betaWithSmallStoppingTolerance =

    2.0999

To initialize the iterations from a qadratic V computed for the linearized dynamics solving A'*P+P*A=-I for P>0, do the following:

opt.sim.flag = 0;
% opt.sim.flag = 1 to use simulation data and =0 with no simulations.
outputs = wrapper(sys,[],roaconstr,opt);
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta =

    2.0985

outputs.RoaEstInfo.info.SimLFG
ans =

     []