Example (E6) 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(1)+x(2)*x(3)^2;-x(2)+x(1)*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'*[0.2849 0.0356 -0.0008; 0.0356 0.2888 -0.0014; -0.0008 -0.0014 0.2835]*x; zV = monomials(x,2:4); Bis.flag = 0; Bis.r1deg = 4;
Generate the options
[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p, Bis); opt.coordoptim.IterStopTol = 1e-4; opt.display.roaest = 1; opt.coordoptim.MaxIters = 30; opt.sim.BetaInit = 2; opt.sim.BetaMax = 2;
Solve the problem
outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations System 1: Num Stable = 100 Num Unstable = 0 Beta for Sims = 2.000 Beta UB = Inf ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 2.000 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 1.401 -------------------Iteration = 1 Beta = 1.612 (Gamma = 0.721) -------------------Iteration = 2 Beta = 1.698 (Gamma = 0.763) -------------------Iteration = 3 Beta = 1.739 (Gamma = 0.787) -------------------Iteration = 4 Beta = 1.767 (Gamma = 0.804) -------------------Iteration = 5 Beta = 1.789 (Gamma = 0.817) -------------------Iteration = 6 Beta = 1.807 (Gamma = 0.827) -------------------Iteration = 7 Beta = 1.821 (Gamma = 0.835) -------------------Iteration = 8 Beta = 1.834 (Gamma = 0.844) -------------------Iteration = 9 Beta = 1.843 (Gamma = 0.850) -------------------Iteration = 10 Beta = 1.851 (Gamma = 0.855) -------------------Iteration = 11 Beta = 1.858 (Gamma = 0.859) -------------------Iteration = 12 Beta = 1.864 (Gamma = 0.863) -------------------Iteration = 13 Beta = 1.869 (Gamma = 0.866) -------------------Iteration = 14 Beta = 1.873 (Gamma = 0.868) -------------------Iteration = 15 Beta = 1.877 (Gamma = 0.870) -------------------Iteration = 16 Beta = 1.880 (Gamma = 0.873) -------------------Iteration = 17 Beta = 1.884 (Gamma = 0.875) -------------------Iteration = 18 Beta = 1.887 (Gamma = 0.877) -------------------Iteration = 19 Beta = 1.890 (Gamma = 0.879) -------------------Iteration = 20 Beta = 1.893 (Gamma = 0.880) -------------------Iteration = 21 Beta = 1.896 (Gamma = 0.882) -------------------Iteration = 22 Beta = 1.898 (Gamma = 0.883) -------------------Iteration = 23 Beta = 1.900 (Gamma = 0.884) -------------------Iteration = 24 Beta = 1.902 (Gamma = 0.886) -------------------Iteration = 25 Beta = 1.904 (Gamma = 0.887) -------------------Iteration = 26 Beta = 1.906 (Gamma = 0.888) -------------------Iteration = 27 Beta = 1.908 (Gamma = 0.890) -------------------Iteration = 28 Beta = 1.910 (Gamma = 0.891) -------------------Iteration = 29 Beta = 1.912 (Gamma = 0.892) -------------------Iteration = 30 Beta = 1.914 (Gamma = 0.893)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 1.9143
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 20
Volume of the ROA estimate
Volume = vol/(4*pi)*3
Volume = 34.6967