% -*-octave-*-
       
clear all;
        
% continuous plant specification: 2 real poles
        
gs = zpk( [], [-1,-2], 1);
        
% convert to continuous state space model for simulation
        
[A,B,C,D] = tf2ss(gs);
        
% design digital controller
        
Ts = 1; % sampling period [s]
gz = c2d( gs, Ts); % convert plant to digital (zoh)
z  = tf('z', Ts); % z operator
T  = 1/z; % dead-beat controller target
kz = T/(gz*(1-T)) % control design!
       
% Control law will be of the form:
% u_k = Kb(1)*e_k + Kb(2)*e_{k-1} ... - Ka(1)*u_{k-1} - Ka(2)*u_{k-2} ...
% copy controller coefficients into matrices for implementation

[b,a] = ss2tf(ss(kz)) % get coefficients of controller
Kb = zeros(1,10); % initialise arrays to zero
Ka = zeros(1,10); % (10 elements is larger than needed)
Kb(1:length(b)) = b./a(1); % populate with controller coefficients
Ka(1:length(a)-1) = a(2:end)./a(1);

% initialise past I/O signals to zero
ez = zeros(size(Ka'));
uz = zeros(size(Kb')); 

% initialise plant and set point
x  = zeros(size(A,1),1);	% zero all states (number from state matrix)
u0  = 0; % zero controller output
ref = 1; % set-point reference

% simulation parameters
dt  = 0.1; % simulation step size
t1  = 10; % simulation end time

% initialise data record (more efficient than dynamic memory reallocation)
data = zeros(t1*dt+1, 3); % number of sample times, [t,u,y]
row = 0; % counter for inserting data records

for t = 0.0 : dt : t1  % loop over time
    row = row + 1;
    y = C*x + D*u0;		% update outputs (continuous)

    % if at sampling time ...
    if (mod(t,Ts) == 0) % update error and control output
       disp('Sampling:'), t
       e0    = ref-y; % measure current error
       ez    = shift(ez,1);	% shift old errors by 1 sample
       ez(1) = e0;		% update register with current error
       u0    = Kb*ez - Ka*uz;	% compute control output!
       uz    = shift(uz,1);	% shift old inputs back 1 sample
       uz(1) = u0; % update register with current output
    end

    % update plant every step
    xdot = A*x + B*u0;	       % update rates using last held input
    x    = x + xdot*dt;	       % Forward Euler integration
    data(row,:) = [t,u0,y];    % save data (continuous)
end

t  = data(:,1);
u  = data(:,2);
y  = data(:,3);

figure(2);
plot(t,y ,'-k', t,u,'xr');
legend({'y(t)','u(k)'})
title('Dead-beat unit step response of continuous plant: G(s) = 1/(s^2+3s+2)')
xlabel('Time /s');
ylabel('Controller and system output');
print -depsc 'deadbeat-continuous.eps'