Difference between revisions of "Lotka Volterra fishing problem (TomDyn/PROPT)"
From mintOC
FelixMueller (Talk | contribs) |
|||
Line 1: | Line 1: | ||
+ | Below you can find the MATLAB file that was used in [[:Category: TomDyn/PROPT | TomDyn/PROPT]] to create the reference solution and its plot for the [[Lotka Volterra fishing problem]]. | ||
+ | |||
+ | |||
<source lang="matlab"> | <source lang="matlab"> | ||
Latest revision as of 18:52, 31 January 2016
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));