% predpreyfbk.m - feedback control of the predator prey systems
% RMM, 13 Sep 06

aminit;					% initialize macros

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

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

% Linearization of the dynamics (from Mma)
xeq = [20.5882; 29.481]; 
App = [0.126275, -0.933333; 0.568, 0];
Bpp = [17.1972; 0];
Cpp = [0 1];

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

% Now set the parameters for our feedback function
predprey_K = Kpp;
predprey_xd = xeq;
predprey_ud = -ref / (Cpp * ((App - Bpp * Kpp) \ Bpp));

% Initial condition and simulation time
x0 = [15, 20];
Tmax = 100;

% Generate a simulation of the system
figure(1); clf; subplot(221);
[t, x] = ode45('predprey', [0 Tmax], x0, [], -r, d, b, k, a, c);
plot(t, x(:,1), '-', t, x(:,2), '--');
xlabel('Time (years)');
ylabel( 'Population', 'Rotation', 90);
amaxis([0 Tmax 0 80]);

lgh = legend('Hare', 'Lynx');
legend(lgh, 'boxoff');
amprint('predprey-sim.eps');

% Save the final point we reached
xfin = x(length(t), :)

% Generate the phase portrait
figure(2); clf; subplot(221); hold on;
amphaseplot('predprey', [0 50 5], [0 50 6], 0.2, [], ...
 0, -r, d, b, k, a, c); 
amphaseplot('predprey', [0 50 5], [60 100 4], 0.13, [], ...
 0, -r, d, b, k, a, c);
amphaseplot('predprey', [60 100 4], [0 50 6], 0.15, [], ...
 0, -r, d, b, k, a, c);
amphaseplot('predprey', [60 100 4], [60 100 4], 0.1, [], ...
 0, -r, d, b, k, a, c);
ah = amphaseplot('predprey', ...
  [0 100 10], [0 100 10], 0, ...
  [ 0.5 10; 100 45; 100 70 ], ...
  Tmax/2, -r, d, b, k, a, c);
set(ah, 'LineWidth', AM_data_linewidth);
xlabel('Hares'); ylabel('Lynxes', 'Rotation', 90);
amaxis([0 100 0 100]); set(gca, 'DataAspectRatio', [1, 1, 1]);

hold on; plot(xfin(1), xfin(2), 'k.', 'MarkerSize', AM_data_markersize*3);
amprint('predprey-pp.eps');
