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