Lotka Volterra fishing problem (switch)

% Copyright 2014 Mathieu Claeys, http://mathclaeys.wordpress.com/

%% Data

% Scaling factors: it is essential to scale variables so that everything % falls within a unit box, otherwise moment problem will be poorly scaled % (Consider moments of the Dirac measure located at 10...) xscale = 1; tscale = 12; %lscale = sqrt(5);

% Relaxation order order = 5; % Choose order >= 1;

%% Create problem definition structure ocpDef.nModes = 2; ocpDef.nStates = 3; ocpDef.nControls = 0; ocpDef.nLifts = 0; ocpDef.scaling.x = [xscale,xscale,xscale]; ocpDef.scaling.t = tscale; ocpDef.scaling.u = 1; ocpDef.scaling.l = 1; ocpDef.dynamics{1} = @(t,x,u,l) [(x(1)-x(1)*x(2)-0.4*x(1))*tscale;

                               (-x(2)+x(1)*x(2)-0.2*x(2))*tscale;
                               ((x(1)-1)^2+(x(2)-1)^2)*tscale]; %[(1-lscale*l(1))*tscale/xscale;(lscale*l(1)-lscale*l(2))*tscale/xscale];

ocpDef.dynamics{2} = @(t,x,u,l) [(x(1)-x(1)*x(2))*tscale;

                               (-x(2)+x(1)*x(2))*tscale;
                               ((x(1)-1)^2+(x(2)-1)^2)*tscale];

ocpDef.runningCost{1} = @(t,x,u,l) 0; ocpDef.runningCost{2} = @(t,x,u,l) 0; ocpDef.initialCost = @(t,x,l) 0; ocpDef.terminalCost = @(t,x,l) x(3); ocpDef.runningConstraints = @(t,x,u,l) [t*(1-t)>=0; % normalized time in [0,1]

                                       x(1) >= 0; % levels are positive
                                       x(2) >= 0;
                                       x(3) >=0;];
                                      % 0==(lscale*l(1))^2-xscale*x(1); %algebraic lifts for square roots
                                       %l(1)>=0;
                                       %0==(lscale*l(2))^2-xscale*x(2);
                                       %l(2)>=0;]; 

ocpDef.initialConstraints = @(t,x,l) [x(1)==0.5;x(2)==0.7;x(3)==0;t==0]; ocpDef.terminalConstraints = @(t,x,l) [t==1]; % NB: to satisfy Putinar's theorem, there should be for each measure a ball % constraint on all variables. Ignore this assumption at your own risks. ocpDef.integralConstraints = {};

%% Construct and solve GloptiPoly problem % (requires the installation of Yalmip and a SDP solver. We recommend Mosek % for speed and SeDuMi for accuracy)

% Create GloptiPoly measure objects with default names (NB: resets all % GloptiPoly states to zero). To impose specific names, construct measures % by hand. measureSystem = switchedMeasureSystem(ocpDef);

% Create GloptiPoly msdp object, modeling a given order moment relaxation P = switchedRelaxation( ocpDef, measureSystem, order );

% Solve GloptiPoly problem %mset('yalmip',true); % default SDP solver of Yalmip will be called. If %commented, Yalmip is not called and GloptiPoly looks for SeDuMi. [status,obj,m] = msol(P); obj

%% Solution extraction npoints.t = 80; npoints.var = 80; % NB: extract moments of one order less, because higher order moments are % very unprecise, since they are not constrained enough [t,x,u,l,d] = extractSolution( ocpDef, measureSystem, order-1, npoints );

figure plot(t,x,'*-'); xlabel('t'); ylabel('x'); legend('x_1','x_2','x_3');

figure plot(t,d(:,1),'o-',t,d(:,2),'+-'); xlabel('t'); ylabel('duty cycle'); legend('Mode 1','Mode 2');

% create/update Bocop directory %toBocop('initLotka',ocpDef,t,x,u,l,d);