Example (E4) from the paper (*)
(*) Topcu, Packard, Seiler, "Local stability analysis using simulations and sum-of-squares programming," Automatica, 2008.
Form the vector field
pvar x1 x2 x3; x=[x1;x2;x3]; f = [-x(2);-x(3);-.915*x(1)+(1-0.915*x(1)^2)*x(2)-x(3)];
p, the basis vector for V and the iteration options (if not equal to the default values) need to be specified before generating the options structures
p = x'*[12.5000 -8.1000 3.0000 -8.1000 20.8000 -8.5000 3.0000 -8.5000 13.4000]*x; zV = monomials(x,2:4); Bis.flag = 0; Bis.r1deg = 4;
Generate the options
[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p, Bis);
Modify the options
opt.coordoptim.IterStopTol = 1e-4; opt.display.roaest = 1;
Solve the problem
outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations System 1: Num Stable = 100 Num Unstable = 0 Beta for Sims = 5.000 Beta UB = Inf ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 5.000 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 4.720 -------------------Iteration = 1 Beta = 5.678 (Gamma = 1.052) -------------------Iteration = 2 Beta = 6.024 (Gamma = 1.075) -------------------Iteration = 3 Beta = 6.174 (Gamma = 1.084) -------------------Iteration = 4 Beta = 6.248 (Gamma = 1.089) -------------------Iteration = 5 Beta = 6.294 (Gamma = 1.092) -------------------Iteration = 6 Beta = 6.318 (Gamma = 1.094) -------------------Iteration = 7 Beta = 6.334 (Gamma = 1.096) -------------------Iteration = 8 Beta = 6.345 (Gamma = 1.097) -------------------Iteration = 9 Beta = 6.353 (Gamma = 1.098) -------------------Iteration = 10 Beta = 6.359 (Gamma = 1.099) -------------------Iteration = 11 Beta = 6.364 (Gamma = 1.099) -------------------Iteration = 12 Beta = 6.367 (Gamma = 1.099) -------------------Iteration = 13 Beta = 6.369 (Gamma = 1.099) -------------------Iteration = 14 Beta = 6.371 (Gamma = 1.100) -------------------Iteration = 15 Beta = 6.373 (Gamma = 1.100) -------------------Iteration = 16 Beta = 6.375 (Gamma = 1.100) -------------------Iteration = 17 Beta = 6.376 (Gamma = 1.100) -------------------Iteration = 18 Beta = 6.377 (Gamma = 1.100) -------------------Iteration = 19 Beta = 6.378 (Gamma = 1.100) -------------------Iteration = 20 Beta = 6.379 (Gamma = 1.100)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 6.3790
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 20
Volume of the ROA estimate
Volume = vol/(4*pi)*3
Volume = 0.6636