Example (E2) 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; x=[x1;x2]; f = [x(2);-2*x(1)-x(2)+x(1)*x(2)^2-x(1)^5+x(1)*x(2)^4+x(2)^5];
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
zV = monomials(x,2:4); Bis.flag = 1; Bis.s2deg = 2; p = x'*[0.9937 0.2283 0.2283 0.8907]*x;
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 = 0 Num Unstable = 1 Beta for Sims = 0.840 Beta UB = 0.840 System 1: Num Stable = 100 Num Unstable = 1 Beta for Sims = 0.840 Beta UB = 0.840 ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 0.840 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 0.547 -------------------Iteration = 1 Beta = 0.735 (Gamma = 0.796) -------------------Iteration = 2 Beta = 0.787 (Gamma = 0.830) -------------------Iteration = 3 Beta = 0.787 (Gamma = 0.830)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 0.7871
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 0.8396
Volume of the ROA estimate
Volume = vol/pi
Volume = 0.9997