% L4_2_predprey.m - MATLAB source code for Lecture 4.2
% RMM, 23 Oct 2012
%
% Required files: predprey.m, predprey_rh.m

% Define global functions for use in the feedback
global predprey_K predprey_xd predprey_ud;

% Set paramater values for the simulation
alpha = 12.5;		% multiplier for hare population
beta = 25;		% multipler for lynx population
gamma = 0.8;		% multiplier for time scale

r = 2 * gamma;			% birth rate of hares
d = 0.7 * gamma;		% death rate for lynxes
b = 0.3 * beta / alpha;		% birth coefficient for lynxes
k = 10 * alpha;			% carrying capacity for hares
a = 8 * gamma * alpha / beta;	% interaction coefficient
c = 4 * alpha;			% saturation coefficient

% Compute the equilibrium point (Note: parameters must be set in predprey.m)
f = inline('predprey(0, x)', 'x');
xeq = fsolve(f, [20, 30])'; He = xeq(1); Le = xeq(2);

% Generate the linearization around the equilibrium point
% App, Bpp computed from symbolic expressions (using Mathematica)
App = [
  -((a*c*k*Le + (c + He)^2*(2*He - k)*r)/((c + He)^2*k)), -((a*He)/(c + He));
  (a*b*c*Le)/(c + He)^2, -d + (a*b*He)/(c + He)
];
Bpp = [He*(1 - He/k); 0];
Cpp = [0 1];

% Assign the desired eigenvalues for the system
Kpp = place(App, Bpp, [-0.1; -0.2]);

% Compute the reference gain
krpp = -1 / (Cpp * ((App - Bpp * Kpp) \ Bpp));

% Now set the parameters for our feedback function
predprey_K = Kpp;
predprey_xd = xeq;

%%
%% Predator prey tracking
%%

figure(3); clf;

% Set initial reference point to equilibrium value
ref = 0; predprey_ud = krpp * ref;
x0 = [20, 30];
T1 = 40;
[t, x] = ode45('predprey', [0 T1], x0, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--'); hold on;
plot([0, T1], [ref + Le, ref+Le], '-k');

% Change the reference point and simulate a bit more
ref = +5; predprey_ud = krpp * ref;
x1 = x(length(t), :); T2 = T1 + 40;
[t, x] = ode45('predprey', [T1 T2], x1, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--'); hold on;
plot([T1, T2], [ref + Le, ref+Le], '-k');

% Change the reference point and simulate a bit more
ref = +10; predprey_ud = krpp * ref;
x2 = x(length(t), :); T3 = T2 + 40;
[t, x] = ode45('predprey', [T2 T3], x2, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--'); hold on;
plot([T2, T3], [ref + Le, ref+Le], '-k');

% Change the reference point and simulate a bit more
ref = 0; predprey_ud = krpp * ref;
x3 = x(length(t), :); T4 = T3 + 40;
[t, x] = ode45('predprey', [T3 T4], x3, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--'); hold on;
plot([T3, T4], [ref + Le, ref+Le], '-k');

% Change the reference point and simulate a bit more
ref = -5; predprey_ud = krpp * ref;
x4 = x(length(t), :); T5 = T4 + 40;
[t, x] = ode45('predprey', [T4 T5], x4, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--'); hold on;
plot([T4, T5], [ref + Le, ref+Le], '-k');

axis([0, 200, 15, 45])
xlabel('Time (years)');
ylabel( 'Population', 'Rotation', 90);

lgh = legend('Hare', 'Lynx');
legend(lgh, 'boxoff');
print -dpng 'predprey-track.png'

