% deg(V) = 6
clear all
pvar x1 x2
x = [x1;x2];
load VDPdataPost
V = VDPdata.Analysis.Deg6.VAfter(end);
p = 0.1*(x1^2+x2^2);

% {x : p(x)<= 1} \subseteq {x : V(x)<=1}?
% Plot the level sets.
[CV,hV] = pcontour(V,1,[-3 3 -3 3]);
set(hV,'color','k','linewidth',2);
hold on;
[cp,hp] = pcontour(p,1,[-3 3 -3 3]);
set(hp,'color','r','linewidth',2);

% Try a scalar multiplier
zS = monomials(x,0:0);
S = polydecvar('s',zS,'vec');
sosconstr{1} = S;
sosconstr{2} = 1-V - S*(1-p);
[info,dopt,sossol] = sosopt(sosconstr,x);
info.feas

% Try a quadratic multiplier
zS = monomials(x,0:2);
S = polydecvar('s',zS,'vec');
sosconstr{1} = S;
sosconstr{2} = 1-V - S*(1-p);
[info,dopt,sossol] = sosopt(sosconstr,x);
info.feas

% Try a quartic multiplier
zS = monomials(x,0:4);
S = polydecvar('s',zS,'vec');
sosconstr{1} = S;
sosconstr{2} = 1-V - S*(1-p);
[info,dopt,sossol] = sosopt(sosconstr,x);
info.feas
S = subs(S,dopt)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = VDPdata.System.f;
Vdot = jacobian(V,x)*f+1e-6*x'*x;
Vdot.maxdeg
V.maxdeg

% question:
% {x : V(x) <= gamma} \subseteq {x : Vdot(x)<=0}?
% {x : p(x)<= 10} \subseteq {x : V(x)<=1}?
% Plot the level sets.
[cVdot,hVdot] = pcontour(Vdot,0,[-3 3 -3 3]);
set(hVdot,'color','r','linewidth',2);
hold on;
[CV,hV] = pcontour(V,0.07,[-3 3 -3 3]);
set(hV,'color','k','linewidth',2);


% Try a quadratic multiplier
% try with gamma = 0.07, 0.01, 0.003
zS = monomials(x,2:2);
S = polydecvar('s',zS,'vec');
sosconstr{1} = S;
sosconstr{2} = -Vdot - S*(0.003-V);
[info,dopt,sossol] = sosopt(sosconstr,x);
info.feas

% Plot {x : V(x) <= 0.003}.
[CV,hV] = pcontour(V,0.003,[-3 3 -3 3]);
set(hV,'color','b','linewidth',2);

% Try a quartic multiplier
zS = monomials(x,2:4);
S = polydecvar('s',zS,'vec');
sosconstr{1} = S;
sosconstr{2} = -Vdot - S*(0.07-V);
[info,dopt,sossol] = sosopt(sosconstr,x);
info.feas



% pcontain example
[gbnds,sout] = pcontain(Vdot, V, monomials(x,1:2));
gamma = gbnds(1);
[CV,hV] = pcontour(V,gamma,[-3 3 -3 3]);
set(hV,'color','m','linewidth',2);