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 = []