%% 034IN "Fundamentals of Automatic Control" - Introduction to MATLAB  Part 6 - LTI Systems
%%
% This is the sixth MATLAB live script of the collection 034IN "Fundamentals of Automatic Control" - Introduction to MATLAB, devoted to introducing the MATLAB environment and tools for solving practical problems related to the topics of the 034IN course, i.e., performance analysis of dynamic systems, design of control law for continuous-time linear dynamical systems, approximate discretization of a continuous-time control law, etc.
% Use this link to go back to the main live script of the collection.
%
% Table of Contents
%  Objectives
% Time Invariant Linear Dynamic Systems
%     LTI-System Object from the State-Space Description
%         A few Examples
%     LTI-System Object from the Transfer Function
%         Example
%         Alternative Way of Using the tf Command
% Simulation of LTI Systems
%     Example: Forced Output Movement of a Continuous-Time System
%     Example: Free Output & State Movements of a Continuous-Time System
%     Example: Output & State Movements of a Continuous-Time System
%     Example: Computation of the Static Gain
% Summary
%         Back to the Index
%         Back to the Previous Part: Visualization & Programming
%%  Objectives
%
%  The aim of this module is to understand how to create and manipulate LTI systems in MATLAB.
%%
%% Time Invariant Linear Dynamic Systems
%
% AN LTI system can be described
% using a state-space description
% using the transfer function matrix for a MIMO system or a single transfer function for a SISO system
% In MATLAB, it is possible to define a linear time-invariant dynamical system as an object of type numLTI model from any of these descriptions (use help numlti for details ). Using these MATLAB objects, it is possible to analyze the properties of the corresponding dynamic system, and also to simulate the evolution of the dynamic system over time with assigned initial conditions and inputs.
% From now on, for simplicity, we only consider SISO systems.
%
%
%% LTI-System Object from the State-Space Description
%
% Starting from the state equation, you can
% define the matrices A, B, C, D of the state equations in the MATLAB workspace;
% define the LTI system by means of the command ss

help ss

%
%% A few Examples
% Given the following state-space description of a continuous-time dynamical system
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& -x_2(t)+3 u(t)\\
% \dot{x}_2(t)&=& -3 x_1(t)+2 x_2(t)\\
% y(t)&=& 4 x_1(t)+2 u(t)
% \end{array}
% \right.
% the corresponding MATLAB object is given by

clear

A = [0 -1;-3 2]; B = [3;0]; C=[4 0]; D=2; % the system matrices
sysC = ss(A, B, C, D) %#ok<NASGU> % the LTI system

%
% Consider a continuous-time dynamical system, described by the model
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& x_2(t)\\
% \dot{x}_2(t)&=&  x_3(t)\\
% \dot{x}_3(t) &=& x_1(t)-2x_2(t)+3x_3(t)+6 u(t)\\
% y(t)&=&  x_1(t)
% \end{array}
% \right.
% then

A = [0 +1 0;0  0  +1; 1  -2 +3]; B = [0; 0; 6]; C=[1  0  0]; D=0; % the system matrices
sysC_2 = ss(A, B, C, D) %#ok<NASGU> % the LTI system

%
%% LTI-System Object from the Transfer Function
%
% Given a transfer function as model for an LTI system, you can create an LTI-model object by using the tf or the zpk command:

help tf
help zpk

%
%% Example
%
% Consider the LTI systems
% G_1(s) = \displaystyle \frac{s+1}{s^2+3s+16} \qquad \qquad G_2(s) = 4\,\displaystyle \frac{(s+1.0)\,(s+0.4)}{s\,(s+8)}
%
% then

numG1 = [ 1 1 ]; denG1 = [ 1 3 16 ];  % coeffs of the polynomials N(s) and D(s)
G1 = tf(numG1, denG1)

Zset = [-1.0  -0.4]; % the zeros
Pset = [0  -8];      % the poles
Kval = 4;            % the gain coeff.
G2 = zpk(Zset, Pset, Kval)

%
%% Alternative Way of Using the tf Command
% Continuous-time and discrete-time ''elementary'' transfer functions can be defined, and used as ''building blocks'' in writing the expressions of the actual transfer functions.
%
% Consider
% G_1(s) = e^{-s}\, \displaystyle \frac{s+1}{s^2+3s+16}

s=tf('s');
Gs = exp(-2*s) * (s+1)/(s^2+3*s+2)

%
%
%% Simulation of LTI Systems
%
% The following MATLAB commands are available, both for continuous-time and discrete-time LTI system:
% step: simulation of the step response;
% initial: simulation of the free state and output movements;
% lsim: simulation of the forced movements.
% dcgain: computation of the static gain (DC gain) of an LTI system
%
%% Example: Forced Output Movement of a Continuous-Time System
% Consider the system
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& -x_1(t)+3 u(t)\\
% \dot{x}_2(t)&=& -3 x_1(t)-2 x_2(t)\\
% y(t)&=& 4 x_1(t)+2 u(t)
% \end{array}
% \right.
% Let's determine the forced state and output trajectories corresponding to the input
% u(t) = \big[3+2\,\sin\left(4 \pi\,t\right) \big] \cdot 1(t)
%

clear

A = [-1 0;-3 -2]; B = [3;0]; C = [4 0]; D = 2; % the system matrices
sysC = ss(A, B, C, D) % the LTI system

Tmax = 10; % let's simulate the system behavior in the time interval [0, Tmax]
t = (0:0.001:Tmax); % a vector containing some time instants
u = 2*sin(2*pi*2*t); % the input values -- Notice: we do not apply the Heaviside function!

[y, ~, x] = lsim(sysC, u, t); % evaluate the forced state & output movement

figure;
subplot(2, 1, 1);
plot(t, x, 'LineWidth', 1.5);
grid on;
xlabel('time t [s]'); ylabel('state vars');
legend('x_{1}(t)', 'x_{2}(t)', 'Location', 'best')

subplot(2,1,2);
plot(t,y, 'LineWidth', 1.5); % plotting the results
xlabel('time t [s]'); ylabel('output');
legend('y(t)', 'Location', 'best')

%
%% Example: Free Output & State Movements of a Continuous-Time System
% Consider the system
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& x_2(t)\\
% \dot{x}_2(t)&=&  x_3(t)\\
% \dot{x}_3(t) &=& -x_1(t)-2x_2(t)-3x_3(t)+6 u(t)\\
% y(t)&=&  x_1(t)
% \end{array}
% \right.
% and let's evaluate the free state and output trajectories, starting from the initial state
% x(0) = \left[ \begin{array}{c}
% 0 \\ -2 \\ +5
% \end{array} \right]

A = [0 +1 0;0  0  +1; -1  -2 -3]; B = [0; 0; 6]; C=[1  0  0]; D=0; % the system matrices
sysC_2 = ss(A, B, C, D); % the LTI system

x0 = [    0
         -2
         +5];

[y, ~, x] = initial(sysC_2, x0, t); % evaluate the forced state & output movement

figure;
subplot(2, 1, 1);
plot(t, x, 'LineWidth', 1.5);
grid on;
xlabel('time t [s]'); ylabel('state vars');
legend('x_{1}(t)', 'x_{2}(t)', 'x_{3}(t)', 'Location', 'best')

subplot(2,1,2);
plot(t,y, 'LineWidth', 1.5); % plotting the results
xlabel('time t [s]'); ylabel('output');
legend('y(t)', 'Location', 'best')

%
%
%% Example: Output & State Movements of a Continuous-Time System
% Consider the system
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& x_2(t)\\
% \dot{x}_2(t)&=&  x_3(t)\\
% \dot{x}_3(t) &=& -x_1(t)-2x_2(t)-3x_3(t)+6 u(t)\\
% y(t)&=&  x_1(t)
% \end{array}
% \right.
% and let's evaluate the state and output trajectories, starting from the initial state
% x(0) = \left[ \begin{array}{c}
% 0 \\ -2 \\ +5
% \end{array} \right]
% when applying the input
% u(t) = \big[3+2\,e^{-0.1 t}\sin\left(4 \pi\,t\right) \big] \cdot 1(t)

clear

A = [0 +1 0;0  0  +1; -1  -2 -3]; 
B = [0; 0; 6];
C=[1  0  0]; D=0; % the system matrices

sysC = ss(A, B, C, D); % the LTI system

x0 = [    0
         -2
         +5];

Tmax = 10; % let's simulate the system behavior in the time interval [0, Tmax]
t = (0:0.001:Tmax); % a vector containing some time instants
u = 3+2 .* exp(-0.1*t).* sin(2*pi*2*t); % the input values -- Notice: we do not apply the Heaviside function!


% free trajs
[yfree, ~, xfree] = initial(sysC, x0, t); % evaluate the forced state & output movement

% forced trajs
[yforced, ~, xforced] = lsim(sysC, u, t); % evaluate the forced state & output movement

x = xfree + xforced;

y = yfree + yforced;

figure;
subplot(2, 1, 1);
plot(t, x, 'LineWidth', 1.5);
grid on;
xlabel('time t [s]'); ylabel('state vars');
legend('x_{1}(t)', 'x_{2}(t)', 'x_{3}(t)', 'Location', 'best')

subplot(2,1,2);
plot(t,y, 'LineWidth', 1.5); % plotting the results
xlabel('time t [s]'); ylabel('output');
legend('y(t)', 'Location', 'best')

%
%% Example: Computation of the Static Gain
% Consider the LTI system
% \left\{
% \begin{array}{lcl}
% \dot{x}_1(t)&=& x_2(t)\\
% \dot{x}_2(t)&=&  x_3(t)\\
% \dot{x}_3(t) &=& -x_1(t)-2x_2(t)-3x_3(t)+6 u(t)\\
% y(t)&=&  x_1(t)
% \end{array}
% \right.
% and let's compute the static gain, by exploiting the MATLAB function dcgain:

clear

A = [0 +1 0;0  0  +1; -1  -2 -3]; 
B = [0; 0; 6];
C=[1  0  0]; D=0; % the system matrices

sysC = ss(A, B, C, D); % the LTI system

sysC_mu = dcgain(sysC)

%
%
% Consider the LTI systems
% G_1(s) = \displaystyle \frac{s+1}{s^2+3s+16} \qquad \qquad G_2(s) = 4\,\displaystyle \frac{(s+1.0)\,(s+0.4)}{s\,(s+8)}
% \qquad \qquad G_3(s) = e^{-s}\, \displaystyle \frac{s+1}{s^2+3s+16}
% and let's compute the static gain:

clear

numG1 = [ 1 1 ]; denG1 = [ 1 3 16 ];  % coeffs of the polynomials N(s) and D(s)
G1 = tf(numG1, denG1)

mu_G1 = dcgain(G1)

Zset = [-1.0  -0.4]; % the zeros
Pset = [0  -8];      % the poles
Kval = 4;            % the gain coeff.
G2 = zpk(Zset, Pset, Kval)

mu_G2 = dcgain(G2)

s=tf('s');
G3 = exp(-2*s) * (s+1)/(s^2+3*s+2)

mu_G3 = dcgain(G3)

%
%% Summary
%
% Using this live script you have learned:
% How to describe an LTI system in MATLAB;
% How to compute in MATLAB the free and the forced state and output movements for an LTI system.
%
%% Back to the Index
% Use this link to go back to the main live script of the collection.
%
%% Back to the Previous Part: Visualization & Programming
% Use this link to go back to the previous live script of this collection.