Example (E7) 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=[-0.42*x(1)-1.05*x(2)-2.3*x(1)^2-0.5*x(1)*x(2)-x(1)^3;1.98*x(1)+x(1)*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 = x'*[29.8951 0.77567/2;0.77567/2 15.727]*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-5; 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.742 -------------------Iteration = 1 Beta = 6.531 (Gamma = 1.259) -------------------Iteration = 2 Beta = 7.091 (Gamma = 1.347) -------------------Iteration = 3 Beta = 7.511 (Gamma = 1.418) -------------------Iteration = 4 Beta = 7.848 (Gamma = 1.476) -------------------Iteration = 5 Beta = 8.121 (Gamma = 1.525) -------------------Iteration = 6 Beta = 8.335 (Gamma = 1.562) -------------------Iteration = 7 Beta = 8.504 (Gamma = 1.591) -------------------Iteration = 8 Beta = 8.639 (Gamma = 1.613) -------------------Iteration = 9 Beta = 8.747 (Gamma = 1.631) -------------------Iteration = 10 Beta = 8.837 (Gamma = 1.645) -------------------Iteration = 11 Beta = 8.912 (Gamma = 1.657) -------------------Iteration = 12 Beta = 8.977 (Gamma = 1.667) -------------------Iteration = 13 Beta = 9.032 (Gamma = 1.676) -------------------Iteration = 14 Beta = 9.078 (Gamma = 1.683) -------------------Iteration = 15 Beta = 9.119 (Gamma = 1.689) -------------------Iteration = 16 Beta = 9.155 (Gamma = 1.695) -------------------Iteration = 17 Beta = 9.185 (Gamma = 1.700) -------------------Iteration = 18 Beta = 9.213 (Gamma = 1.704) -------------------Iteration = 19 Beta = 9.237 (Gamma = 1.707) -------------------Iteration = 20 Beta = 9.257 (Gamma = 1.710)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 9.2570
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 20
Volume of the ROA estimate
Volume = vol/pi
Volume = 1.4737