Robust ROA calculations

dynamics:

x1dot = x2;

x2dot = -x2-2*x1+2*x1^3 + delta*(-x1+x1^3);

with delta \in [-1,1]

This example was also used in Topcu and Packard, IEEE TAC, 2009 (in the special issue on positive polynomials in controls (example 1 in the paper)

% Form the vector field
pvar x1 x2;
x = [x1;x2];
x1dot = x2;
x2dot = -x2-2*x1+2*x1^3;

Nominal system

f = [x1dot; x2dot];

Introduce an uncertain parameter

pvar d1

Specify its range

ini_cell = [-1 1];

Form the uncertain vector field

f = f + d1*[0; -x1+x1^3];


% Get the vertex system
[roaconstr,opt,sys] = GetRoaOpts(f, x);
[fNOM,fVER] = getf(sys,ini_cell);

% Generate the options, etc.
zV = monomials(x,2:4);
Bis.flag = 0;
Bis.r1deg = 4;


[roaconstr,opt,sys] = GetRoaOpts(fVER, x, zV, [], Bis);
sys.fWithDel = [];

opt.sim.NumConvTraj = 40;
opt.display.roaest = 1;

Run the computations

outputs = wrapper(sys,[],roaconstr,opt);
------------------Beginning simulations
System 1: Num Stable = 0 	Num Unstable = 1 	Beta for Sims = 3.289 	Beta UB = 3.289 
System 1: Num Stable = 0 	Num Unstable = 2 	Beta for Sims = 1.390 	Beta UB = 1.390 
System 1: Num Stable = 2 	Num Unstable = 3 	Beta for Sims = 1.306 	Beta UB = 1.306 
System 1: Num Stable = 4 	Num Unstable = 4 	Beta for Sims = 0.913 	Beta UB = 0.913 
System 1: Num Stable = 6 	Num Unstable = 5 	Beta for Sims = 0.861 	Beta UB = 0.861 
System 1: Num Stable = 12 	Num Unstable = 6 	Beta for Sims = 0.818 	Beta UB = 0.842 
System 1: Num Stable = 18 	Num Unstable = 7 	Beta for Sims = 0.777 	Beta UB = 0.808 
System 2: Num Stable = 1 	Num Unstable = 1 	Beta for Sims = 1.476 	Beta UB = 0.808 
System 2: Num Stable = 3 	Num Unstable = 2 	Beta for Sims = 1.402 	Beta UB = 0.808 
System 2: Num Stable = 6 	Num Unstable = 3 	Beta for Sims = 1.114 	Beta UB = 0.808 
System 2: Num Stable = 6 	Num Unstable = 4 	Beta for Sims = 1.058 	Beta UB = 0.808 
System 2: Num Stable = 8 	Num Unstable = 5 	Beta for Sims = 1.000 	Beta UB = 0.808 
System 2: Num Stable = 10 	Num Unstable = 6 	Beta for Sims = 0.929 	Beta UB = 0.808 
System 2: Num Stable = 11 	Num Unstable = 7 	Beta for Sims = 0.882 	Beta UB = 0.808 
------------------End of simulations
------------------Begin search for feasible V
Try = 1 	 Beta for Vfeas = 0.882
Try = 2 	 Beta for Vfeas = 0.838
------------------Found feasible V
Initial V (from the cvx outer bnd) gives Beta = 0.173
-------------------Iteration = 1 
Beta = 0.567 (Gamma = 0.535) 
-------------------Iteration = 2 
Beta = 0.665 (Gamma = 0.604) 
-------------------Iteration = 3 
Beta = 0.716 (Gamma = 0.640) 
-------------------Iteration = 4 
Beta = 0.739 (Gamma = 0.656) 

Extract the solution

[V,beta,gamma,p,multip,betaUpper] = extractSol(outputs);
beta
beta =

    0.7388

Upper bound on beta

betaUpper
betaUpper =

    0.8822

Plot the results

[Cp4,hp4] = pcontour(p,beta,[-2 2 -2 2],'k'); hold on;
set(hp4,'linewidth',2);
[CV4,hV4] = pcontour(V,gamma,[-2 2 -2 2],'b');
set(hV4,'linewidth',2);
set(gca,'xlim',[-1.5 1.5],'ylim',[-1.5 1.5]);

traj = outputs.RoaEstInfo.info.SimLFG.sim.Trajectories(1).unstab(end).state;
pval = peval(traj,p.coef,p.deg);
[aux,ind] = min(pval);
plot(traj(1,ind),traj(2,ind),'r*','markersize',8);
grid on;