% Drug administration example
% RMM, 13 Sep 06

aminit;					% initialize macros

% Parameter settings for the model (need comments)
k0=0.1;
k1=0.1;
k2=0.5;
b0=1.5;

% Compartment model definition
Ada = [-k0-k1 k1; k2 -k2];
Bda = [b0;0];
Cda = [0 1];
sysol = ss(Ada, Bda, Cda, 0);

% Determine the desired eigenvalue locations based on desired dynamics
zeta = 0.9;
omega0 = 2.2/(10);
lamr = -zeta*omega0;
lami = omega0*sqrt(zeta^2 - 1);

% Determine where to set the gains
Kda = place(Ada, Bda, [lamr+lami, lamr-lami]);

% Close the loop
kr = - 1/(Cda * ((Ada - Bda * Kda) \ Bda));
syscl = ss(Ada - Bda*Kda, Bda * kr, Cda, 0);

% Simulation parameters
Tper = 0.01;				% simulation period
Tdur = 1.04;				% pulse duration
Tmax = 50;				% simulation length

% Set up a pulse train for the open loop system
pulse = (Tper:Tper:10) * 0; pulse(1:Tdur/Tper) = 1;
pills = [pulse pulse pulse pulse pulse 0] / 1.5;
tvec = 0:Tper:Tmax;

% Simulate open loop and closed loop
[yol, tolv, xol] = lsim(sysol, pills, tvec);
[ycl, tclv, xcl] = step(syscl, Tmax);

figure(1); clf; subplot(311);
plot(tclv, ycl, '-', tolv, yol, '--', 'Linewidth', AM_data_linewidth);
amaxis([0 50 0 1.5]);
xlabel('Time (min)'); ylabel('Concentration, C2');
lgh = legend('State feedback', 'Pulses', 'Location', 'SouthEast');
legend(lgh, 'boxoff');

subplot(312);
plot(tclv, kr - xcl*Kda', '-', tvec, pills, '--', ...
  'Linewidth', AM_data_linewidth);
amaxis([0 50 0 0.8]);
xlabel('Time (min)'); ylabel('Input dosage');

amprint('compartment-statefbk.eps');
