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