Robust ROA calculations for the Van der Pol dynamics

Contents

This file demonstrates the problem setup for robust ROA computation with parametric uncertainty. For more details on all the options see NACodeExplanation.pdf.

Affine uncertainty

Form the vector field

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

Nominal system

f = [x1dot; x2dot];

Introduce an uncertain parameter

pvar d1

Specify its range

ini_cell = [-0.1 0.1];

Form the uncertain vector field

f = f + d1*[-x2;0];

Set some options

[roaconstr,opt,sys] = GetRoaOpts(f, x);
opt.sim.NumConvTraj = 40;
Solve the problem for an estimate of the robust ROA
outputs = wrapper(sys,ini_cell,roaconstr,opt);

Extract the solution

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

Certified beta

beta
beta =

    1.3030

Upper bound on beta

betaUpper
betaUpper =

    2.3544

The above code implement the two-step procedure for robust ROA calculations:

- Compute an "optimal" V for the center of the cell - Assess V from the initial step for the vertex systems

We can try to compute an optimal V directly taking all the vertices in to account in constructing V

[fNOM,fVER] = getf(sys,ini_cell);

Re-generate the options, etc.

[roaconstr,opt,sys] = GetRoaOpts(fVER, x);
sys.fWithDel = [];

Run the computations

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

Extract the solution

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

    1.3256

Non-affine uncertainty

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

Introduce the uncertain parameters

pvar d1 d2

Specify its value ini_cell = [d1min d1max d2min d2max ... dmmin dmmax]

ini_cell = [-.2 .2 -.1 .2];

Form the uncertain vector field

f = f + d1*[-x2;0] + d2*[0;-x1] + d1^2*[-x1+x2;0]
[roaconstr,opt,sys] = GetRoaOpts(f, x);
opt.sim.NumConvTraj = 40;
 
f = 
  [ -d1^2*x1 + d1^2*x2 - d1*x2 - x2 ]
  [       x1^2*x2 - d2*x1 + x1 - x2 ]
 

Solve

outputs = wrapper(sys,ini_cell,roaconstr,opt);

Extract the solution

[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs)
 
V = 
  0.47208*x1^2 - 0.19451*x1*x2 
  + 0.43122*x2^2
 

beta =

    0.7321


gamma =

    0.4034

 
p = 
  x1^2 + x2^2
 

multip = 

    r1: [8x1 polynomial]
    r2: [1x1 polynomial]


betaUpper =

    2.1455