Lotka Volterra fishing problem (switch)

Below you can find the MATLAB file that was used in switch to create the reference solution and its plot for the Lotka Volterra fishing problem.

% 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);