Example (E3) 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);x(3);-4*x(1)-3*x(2)-3*x(3)+x(1)^2*x(2)+x(1)^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.6352 0.2056 0.0745 0.2056 0.2763 0.0778 0.0745 0.0778 0.0402]*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-3; opt.display.roaest = 1; opt.sim.BetaMax = 1; opt.sim.BetaInit = 1;
Solve the problem
outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations System 1: Num Stable = 17 Num Unstable = 1 Beta for Sims = 0.950 Beta UB = 0.968 System 1: Num Stable = 100 Num Unstable = 1 Beta for Sims = 0.950 Beta UB = 0.968 System 1: Num Stable = 128 Num Unstable = 2 Beta for Sims = 0.902 Beta UB = 0.923 ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 0.902 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 0.447 -------------------Iteration = 1 Beta = 0.568 (Gamma = 0.668) -------------------Iteration = 2 Beta = 0.669 (Gamma = 0.760) -------------------Iteration = 3 Beta = 0.712 (Gamma = 0.810) -------------------Iteration = 4 Beta = 0.736 (Gamma = 0.842) -------------------Iteration = 5 Beta = 0.754 (Gamma = 0.867) -------------------Iteration = 6 Beta = 0.771 (Gamma = 0.891) -------------------Iteration = 7 Beta = 0.781 (Gamma = 0.905) -------------------Iteration = 8 Beta = 0.791 (Gamma = 0.919) -------------------Iteration = 9 Beta = 0.796 (Gamma = 0.926) -------------------Iteration = 10 Beta = 0.800 (Gamma = 0.932) -------------------Iteration = 11 Beta = 0.804 (Gamma = 0.938) -------------------Iteration = 12 Beta = 0.807 (Gamma = 0.943) -------------------Iteration = 13 Beta = 0.810 (Gamma = 0.947) -------------------Iteration = 14 Beta = 0.812 (Gamma = 0.951) -------------------Iteration = 15 Beta = 0.815 (Gamma = 0.954) -------------------Iteration = 16 Beta = 0.816 (Gamma = 0.956) -------------------Iteration = 17 Beta = 0.818 (Gamma = 0.959) -------------------Iteration = 18 Beta = 0.820 (Gamma = 0.962) -------------------Iteration = 19 Beta = 0.821 (Gamma = 0.964) -------------------Iteration = 20 Beta = 0.823 (Gamma = 0.966)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta = 0.8229
betaUpper [vol,volstd] = pvolume(V,gamma,[],10000);
betaUpper = 0.9229
Volume of the ROA estimate
Volume = vol/(4*pi)*3
%
Volume = 21.4432