Example (E1) 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)-3*x(2)+x(1)^2*x(2)];
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 = 0.15515*x1^2 + 0.056512*x1*x2 + 0.055778*x2^2; 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 = 0 Num Unstable = 1 Beta for Sims = 0.923 Beta UB = 0.923 System 1: Num Stable = 9 Num Unstable = 2 Beta for Sims = 0.877 Beta UB = 0.922 System 1: Num Stable = 100 Num Unstable = 2 Beta for Sims = 0.877 Beta UB = 0.922 ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 0.877 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 0.632 -------------------Iteration = 1 Beta = 0.789 (Gamma = 0.770) -------------------Iteration = 2 Beta = 0.805 (Gamma = 0.783) -------------------Iteration = 3 Beta = 0.813 (Gamma = 0.790) -------------------Iteration = 4 Beta = 0.818 (Gamma = 0.795) -------------------Iteration = 5 Beta = 0.822 (Gamma = 0.799) -------------------Iteration = 6 Beta = 0.826 (Gamma = 0.802) -------------------Iteration = 7 Beta = 0.828 (Gamma = 0.804) -------------------Iteration = 8 Beta = 0.831 (Gamma = 0.806) -------------------Iteration = 9 Beta = 0.833 (Gamma = 0.808) -------------------Iteration = 10 Beta = 0.835 (Gamma = 0.809) -------------------Iteration = 11 Beta = 0.837 (Gamma = 0.811) -------------------Iteration = 12 Beta = 0.838 (Gamma = 0.812) -------------------Iteration = 13 Beta = 0.840 (Gamma = 0.813) -------------------Iteration = 14 Beta = 0.841 (Gamma = 0.814) -------------------Iteration = 15 Beta = 0.842 (Gamma = 0.814) -------------------Iteration = 16 Beta = 0.842 (Gamma = 0.815) -------------------Iteration = 17 Beta = 0.843 (Gamma = 0.816) -------------------Iteration = 18 Beta = 0.844 (Gamma = 0.816) -------------------Iteration = 19 Beta = 0.844 (Gamma = 0.817) -------------------Iteration = 20 Beta = 0.845 (Gamma = 0.817)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 0.8447
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 0.9224
Volume of the ROA estimate
Volume = vol/pi
Volume = 14.9526