Example (E5) 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)+2*x(2)*x(3);x(3);-0.5*x(1)-2*x(2)-x(3)]; V = .36158496756578887372799283958143e-2*x2*x1*x3^2+.29783211619555115891485746361835e-2*x2^2*x1*x3+.18439685630571241996894822265806e-2*x1^2*x3*x2+.13044603999184172971795243493243e-1*x1*x3*x2+.30007289365477082998844387631829e-2*x3^3-.22810294738224829698229100890553e-1*x2^3+.14217732978748812961905788581153e-2*x1^3+.51212695564921775109857268792288e-2*x3^4+.26846700085802117677246771865053e-1*x2^4+.18885654998064717995930557086088e-3*x1^4+.70958863760435608253437047048148e-2*x2*x1^2+.22471172198359998990629010870964e-2*x1^2*x3^2+.66291488272207317773332665473139e-3*x1^3*x3+.13372522651752905324027401449332e-2*x1^2*x2^2+.95760444516384065160595927127055e-3*x1^3*x2+.10413330077863903122092369360588e-1*x2*x3^2-.12683893481358567384339869741544e-1*x3*x2^2+.10037589955128897301826355035902e-1*x1*x3^2+.28192806154596763983159780274800e-3*x3*x1^2-.27818432359580794547694733274931e-1*x1*x2^2+.16320784603889616447554694939541e-2*x3^3*x1+.14684470527571642009916226738654e-1*x2^2*x3^2+.50256266161596305749180141420360e-2*x2^3*x3+.11622803374754910947908645288631e-2*x2^3*x1+.51179903699664901874634052996679e-2*x3^3*x2+.42353426634689856060385134678601e-1*x2*x1+.23889147105175959137072900566088e-1*x3*x1+.28485200286470550330566787785647e-1*x1^2+.33718235450379904118271641834909e-1*x2^2+.25730410339349391294856467489531e-1*x3*x2+.16987279265235564953466555985293e-1*x3^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'*[0.65707 0.4468 0.36154 0.4468 1.6409 0.37826 0.36154 0.37826 0.83046]*x; zV = monomials(x,2:4); Bis.flag = 1; Bis.s2deg = 2;
Generate the options
[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, p, Bis);
Modify the options
opt.coordoptim.IterStopTol = 1e-3; 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 = 2.695 -------------------Iteration = 1 Beta = 3.125 (Gamma = 0.581) -------------------Iteration = 2 Beta = 3.516 (Gamma = 0.659) -------------------Iteration = 3 Beta = 3.867 (Gamma = 0.737) -------------------Iteration = 4 Beta = 4.258 (Gamma = 0.825) -------------------Iteration = 5 Beta = 4.609 (Gamma = 0.913) -------------------Iteration = 6 Beta = 5.000 (Gamma = 1.006) -------------------Iteration = 7 Beta = 5.352 (Gamma = 1.099) -------------------Iteration = 8 Beta = 5.703 (Gamma = 1.191) -------------------Iteration = 9 Beta = 6.016 (Gamma = 1.270) -------------------Iteration = 10 Beta = 6.328 (Gamma = 1.348) -------------------Iteration = 11 Beta = 6.641 (Gamma = 1.431) -------------------Iteration = 12 Beta = 6.953 (Gamma = 1.514) -------------------Iteration = 13 Beta = 7.422 (Gamma = 1.655) -------------------Iteration = 14 Beta = 7.773 (Gamma = 1.768) -------------------Iteration = 15 Beta = 8.047 (Gamma = 1.855) -------------------Iteration = 16 Beta = 8.750 (Gamma = 2.119) -------------------Iteration = 17 Beta = 8.750 (Gamma = 2.119)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 7.4219
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 20
Volume of the ROA estimate
Volume = vol/(4*pi)*3
Volume = 62.1598