Upper bound demonstrations
Contents
VDP with deg(V) = 2
Form the vector field
pvar x1 x2; x = [x1;x2]; x1dot = -x2; x2dot = x1+(x1^2-1)*x2; f = [x1dot; x2dot];
Get the default values of options to run the ROA code.
zV = monomials(x,2:2); Bis.r1deg = 2;
Now, call GetRoaOpts to get the correponding opts, roaconstr, etc.
[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, [],Bis); opt.sim.NumConvTraj = 200; opt.sim.dispEveryNth = 40; opt.display.roaest = 1; opt.coordoptim.IterStopTol = 1e-4;
Call the wrapper which in turn calls RoaEst.m
outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations System 1: Num Stable = 0 Num Unstable = 1 Beta for Sims = 2.348 Beta UB = 2.348 System 1: Num Stable = 40 Num Unstable = 1 Beta for Sims = 2.348 Beta UB = 2.348 System 1: Num Stable = 43 Num Unstable = 2 Beta for Sims = 2.230 Beta UB = 2.348 System 1: Num Stable = 80 Num Unstable = 2 Beta for Sims = 2.230 Beta UB = 2.348 System 1: Num Stable = 120 Num Unstable = 2 Beta for Sims = 2.230 Beta UB = 2.348 System 1: Num Stable = 160 Num Unstable = 2 Beta for Sims = 2.230 Beta UB = 2.348 System 1: Num Stable = 200 Num Unstable = 2 Beta for Sims = 2.230 Beta UB = 2.348 ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 2.230 Try = 2 Beta for Vfeas = 2.119 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 1.496 -------------------Iteration = 1 Beta = 1.513 (Gamma = 0.746) -------------------Iteration = 2 Beta = 1.516 (Gamma = 0.747) -------------------Iteration = 3 Beta = 1.517 (Gamma = 0.747) -------------------Iteration = 4 Beta = 1.517 (Gamma = 0.747)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
Upper bounds from divergent trajectories
betaUpperDivergent = outputs.RoaEstInfo.info.SimLFG.sim.BetaUB; betaUpperDivergent
betaUpperDivergent = 2.3478
Upper bound from infeasibility of the relaxation
if betaUpper < outputs.RoaEstInfo.info.SimLFG.sim.BetaUB betaUpperInfeas = betaUpper; betaUpperInfeas else display('No upper bound from infeasibility'); end
betaUpperInfeas = 2.2304
Certified beta
betaCertified = beta; betaCertified
betaCertified = 1.5168
VDP with deg(V) = 6
Form the vector field
pvar x1 x2; x = [x1;x2]; x1dot = -x2; x2dot = x1+(x1^2-1)*x2; f = [x1dot; x2dot];
Get the default values of options to run the ROA code.
zV = monomials(x,2:6); Bis.r1deg = 4;
Now, call GetRoaOpts to get the correponding opts, roaconstr, etc.
[roaconstr,opt,sys] = GetRoaOpts(f, x, zV, [],Bis); opt.sim.NumConvTraj = 200; opt.sim.dispEveryNth = 40; opt.display.roaest = 1; opt.coordoptim.IterStopTol = 1e-4;
Call the wrapper which in turn calls RoaEst.m
outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations System 1: Num Stable = 0 Num Unstable = 1 Beta for Sims = 2.352 Beta UB = 2.352 System 1: Num Stable = 7 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 System 1: Num Stable = 40 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 System 1: Num Stable = 80 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 System 1: Num Stable = 120 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 System 1: Num Stable = 160 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 System 1: Num Stable = 200 Num Unstable = 2 Beta for Sims = 2.235 Beta UB = 2.347 ------------------End of simulations ------------------Begin search for feasible V Try = 1 Beta for Vfeas = 2.235 ------------------Found feasible V Initial V (from the cvx outer bnd) gives Beta = 0.583 -------------------Iteration = 1 Beta = 1.436 (Gamma = 1.387) -------------------Iteration = 2 Beta = 1.731 (Gamma = 1.618) -------------------Iteration = 3 Beta = 1.886 (Gamma = 1.770) -------------------Iteration = 4 Beta = 1.981 (Gamma = 1.874) -------------------Iteration = 5 Beta = 2.049 (Gamma = 1.948) -------------------Iteration = 6 Beta = 2.098 (Gamma = 2.001) -------------------Iteration = 7 Beta = 2.134 (Gamma = 2.041) -------------------Iteration = 8 Beta = 2.164 (Gamma = 2.072) -------------------Iteration = 9 Beta = 2.188 (Gamma = 2.098) -------------------Iteration = 10 Beta = 2.209 (Gamma = 2.119) -------------------Iteration = 11 Beta = 2.228 (Gamma = 2.137) -------------------Iteration = 12 Beta = 2.244 (Gamma = 2.153) -------------------Iteration = 13 Beta = 2.259 (Gamma = 2.166) -------------------Iteration = 14 Beta = 2.272 (Gamma = 2.178) -------------------Iteration = 15 Beta = 2.283 (Gamma = 2.188) -------------------Iteration = 16 Beta = 2.293 (Gamma = 2.197) -------------------Iteration = 17 Beta = 2.301 (Gamma = 2.204) -------------------Iteration = 18 Beta = 2.308 (Gamma = 2.209) -------------------Iteration = 19 Beta = 2.313 (Gamma = 2.213) -------------------Iteration = 20 Beta = 2.317 (Gamma = 2.217)
Extract the solution
[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
Upper bounds from divergent trajectories
betaUpperDivergent = outputs.RoaEstInfo.info.SimLFG.sim.BetaUB; betaUpperDivergent
betaUpperDivergent = 2.3470
Upper bound from infeasibility of the relaxation
if betaUpper < outputs.RoaEstInfo.info.SimLFG.sim.BetaUB betaUpperInfeas = betaUpper; betaUpperInfeas else display('No upper bound from infeasibility'); end
No upper bound from infeasibility
Certified beta
betaCertified = beta; betaCertified
betaCertified = 2.3173