% r_hw1
% Workshop 1 at MPC Workshop in Trondheim (January 2015)
% function files used: examp3tnk.m (Problem 3)
%                      linodepar.m (Problme 4)
%
% solution to homework 1: MPC - Fall, 2013
% 
% b.w. bequette
% 4 Sept 2013
%
%  consider the state space model xdot = a*x + b*u
%                                    y = c*x + d*u
%  the matrices for the 3-tank example are:
  a = [-1 0 0;1 -1 0;0 1 -1]
  b = [1;0;0]
  c = [0 0 1]
  d = 0
% ---------- Problem 1 ---------------------------------------
% find the eigenvalues of the A matrix
%
  eig(a)
%
% ---------- Problem 2 ---------------------------------------
% define the state space "object"
  lintank = ss(a,b,c,d)
% use the default plot function (starts at time = 0)
  step(lintank)
% (notice that the time should really be minutes, not seconds)
% ---------- Problem 3 ---------------------------------------
% integrate the state space model from t = -5 to t = 0 minutes
%    assuming u = 0, followed by an integration from 
%    t = 0 to t = 20 minutes, with u = 1
%    in this simple solution, the step change is embedded
%    in the function file
  tbeg = -5;       % initial simulation time, minutes
  tend = 20;       % final simulation time, minutes
  x0   = [0;0;0];  % initial conditions, column vector
%
  [t,x] = ode45('examp3tnk',[tbeg tend],x0);
%   returns a time vector and an array of states
%   the states are oriented "column-wise"
%
  figure(2)
  plot(t,x(:,3));  % plots the third state
  xlabel('time, min')
  ylabel('x3, meters')
  title('step change in flowrate, embedded ode45 solution')
% ---------- Problem 4 ----------------------------------------
% same step change as problem 3, but pass the parameters and input
%   into the function file
% 
% set up a for loop and integrate over 0.5 minute time steps
%
  delt = 0.5;  % sample time = 0.5 minutes
%   start the simulation at t = -5
  tbeg = -5;
  tend = 20;
  x0   = zeros(3,1);
  time = tbeg:delt:tend;   % generate the time vector
%      for-loop to make an input step change at t = 0 minutes
  for k = 1:length(time)-1;
        if time(k) < 0;
          uvec(k) = 0;
      else
          uvec(k) = 1;
        end
  end
  xdis(:,1) = x0;  % state initial condition
  for k = 1:length(time)-1;
%          generates tdummy (time vector) and xdummy(states array) for time points from
%          time(k) to time(k+1)
      [tdummy,xdummy] = ode45('linodepar',[time(k) time(k+1)],xdis(:,k),[],a,b,uvec(k));
      ndum = length(tdummy);         % finds length of tdummy vector
      xdis(:,k+1) = xdummy(ndum,:)'; % xdummy is oriented with states in columns
  end
      uvec(k+1) = uvec(k);  % need length of uvec to be same as time vector
%
  figure(5)
  plot(time,xdis(3,:))
  xlabel('time, min')
  ylabel('x3, m')
  title('ode45 solution, parameters passed through argument list')
%     use the subplot command to plot the third state (top) and
%     input (bottom)
  figure(7)
  subplot(2,1,1)
  plot(time,xdis(3,:))
  xlabel('t')
  ylabel('x3, m')
  title('using the subplot command')
  subplot(2,1,2)
  plot(time,uvec)
  xlabel('t, min')
  ylabel('u')
%
%  notice that a discontinuous step looks better
%
  [tt,uu] = stairs(time,uvec);  % uses zero-order hold
%
  figure(8)
  subplot(2,1,1)
  plot(time,xdis(3,:))
  axis([-5 20 -0.1 1.1])   % axis limits for the state 3 plot
  xlabel('t, min')
  ylabel('x3, meters')
  subplot(2,1,2)
  plot(tt,uu)
  axis([-5 20 -0.1 1.1])   % axis limits for the input plot
  xlabel('t, min')
  ylabel('u')
% ------------- Problem 5 -----------------------------------------
% create the continuous time transfer function and calculate the poles
%
  [num,den] = ss2tf(a,b,c,d)
%
  roots(den)      % find the poles (inverse minutes)