Lotka Volterra fishing problem (TomDyn/PROPT)

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

%% Lotka-Volterra fishing example (relaxed version)
%
% Implemented by: Maximilian von Wolff
%
%
%% Problem Description
%
% As on https://modest.math.uni-magdeburg.de/mintoc/index.php/Lotka_Volterra_fishing_problem
%
% \min\limits_{x,u}      x_3(t_f)
% 
%  subject to:
% 
%   dot(x)_1(t)   =      x_1(t) - x_1(t)*x_2(t) - c_0*x_1(t)*u(t),
%   dot(x)_2(t)   =      -x_2(t) + x_1(t)*x_2(t) - c_1*x_2(t)*u(t),
%   dot(x)_3(t)   =      (x_1(t) - 1)^2 + (x_2(t) - 1)^2,
% 
% initial conditions:
%
% x(0)=[0.5, 0.7, 0]'
% 
% control constraints:
%
% u(t) \in \{0, 1\}              \forall t \in [t_0, t_f]
%
% Parameters
%
% [t_0, t_f]    =   [0, 12]
% [c_0, c_1]    =   [0.4, 0.2]
%
%
% 
%
%
%% Problem setup
%
% Setup parameters:
% 
n = 400;
%
 
toms t
p = tomPhase('p', t, 0, 12, n);
setPhase(p);
tomStates x1 x2 x3
tomControls u    %or use integer u and set problem type as MINLP and solver as KNITRO as in annotation below
 
%initial states
xi = [0.5, 0.7, 0];
 
 
%initial guess
x0 = {collocate(u == 0)
    icollocate({x1 == xi(1)
                x2 == xi(2)
                x3 == xi(3)
                })};
 
%Box constraints
cbox = {(0 <= collocate(u) <= 1)};
 
%Boundary constraints
cbnd = {initial({   x1 == xi(1); 
                    x2 == xi(2);
                    x3 == xi(3);})};
 
%Parameters for ODE's
c = [0.4, 0.2]';
c0 = c(1);
c1 = c(2);
 
%ODE's
 
dx1 = x1 - x1.*x2 -c0*x1.*u;
dx2 = -x2 + x1.*x2 - c1*x2.*u;
dx3 = (x1 - 1).^2 + (x2 - 1).^2;
 
 
ceq = collocate({
        dot(x1) == dx1
        dot(x2) == dx2
        dot(x3) == dx3});
 
 
%Objective
 
objective = final(x3);
 
%% Solve the problem
    options = struct;
    options.name = 'Lotka-Volterra';
    %options.type = 'minlp';
    %options.solver = 'knitro';
    solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);
 
    xopt1 = subs(x1,solution);
    xopt2 = subs(x2,solution);
    xopt3 = subs(x3,solution);
    uopt  = subs(u,solution);
 
 
%Plotting solution
 
figure(1)
subplot(3,1,1);
ezplot([x1; x2]); legend('Prey','Predator');
xlabel('time');
ylabel('Biomass');
title('Biomasses');
 
subplot(3,1,2);
ezplot(u); legend('u');
xlabel('time');
ylabel('Control');
title('Fishing control');
 
subplot(3,1,3);
ezplot(x3); legend('Obj. Value');
xlabel('time');
title('Objective');
 
disp('Objective Function Value');
disp(final(xopt3));