Hands On: State Space Descriptions - Solution Proposals

Table of Contents

Introduction

This live script illustrates some models of dynamical systems (linear or nonlinear, time-invariant or time-variant, continuous or discrete time).

Inspired by the examples illustrated in the live script StateSpaceDescriptionExamples.m, you will able to create a Matlab live script to simulate the behavior of a continuous time or discrete time dynamical system.

What you will learn:

Given a state space description of a dynamical system, how to build a MATLAB model for such a system;
Tune the model, setting the model parameters (or the simulation parameters) using a Matlab script;
How to run a MATLAB model and retrieve the results.

A Linear Continuous Time System: a Mechanical System

The dynamics of the system is given by the state-space model derived from conservation balances (of mass, energy, etc.). Applying Newton's second law, we obtain the following state-space description:

where the state variables are

This state-space description may be easily implemented in Matlab.

Proposal: implement a MATLAB LTI model and a Matlab live script to analyse the evolution of the above described mechanical system. Use the following parameters to configure the system model:

point masses [kg]: , and ;
viscous damping coefficients [N s/m]: ;
spring constants [N/m] : ;
initial conditions:
external forces (inputs):
simulation time interval [s]:

A Solution Proposal

System Parameters

clear variables
 
% ---------- the system parameters ----------
% point masses 
M1 = 1.0; % [kg]
M2 = 2.5; % [kg]
 
% viscous damping coefficients
B  = 3.6e-3; % [N sec/m]
B1 = 5.0e-2; % [N sec/m]
B2 = 5.0e-2; % [N sec/m]
 
% spring constants
K  = 9.1e-2; % [N/m]
K1 = 1.0e-1; % [N/m]
K2 = 1.0e-1; % [N/m]
% ---------- the system parameters ----------
 
% ---------- the initial state --------------
% x1 <--> position of the mass M1
% x2 <--> speed of mass M1
% x3 <--> position of the mass M2
% x4 <--> speed of the mass M2
x0IC = [ 1.0 ;
           0 ;
        -0.5 ;
           0  ];
% ---------- the initial state --------------
 
% ----------- simulation time interval -------------------------------
% Let's create a "dense" array of time instants. We will determine 
% both the free state movement (the zero-input movement) and the forced
% state movement (the zero-state movement) in these time instants.
 
N = 20000; % how many time instants?
Tstart = 0; % the starting time instant of the simulation
Tstop = 30; % the simulation must end when the time reaches this value
 
tSET = linspace(Tstart, Tstop, N); 
% The array containing the time instants
% at which the values of free and forced
% movements will be determined.
% ----------- simulation time interval -------------------------------

Inputs Generation

Now let's generate an array of values for each external input:

% ----------------- generation of the input sequences --------------------
% the first input: a step signal, with assigned amplitude,
% applied at the time instant t=0
 
U1 = 2; % the amplitude of the 1st input
u1SEQ = U1 * ones(size(tSET));
% the array containing the samples of the input u1(t):
% one sample for each time instant in tSET
 
% the second input: a delayed step input, with assigned amplitude,
% applied starting from the time instant t = 10 sec
U2 = -2;     % the amplitude of the 2nd input
tDELAY = 10; % the starting instant for the delayed step input
 
% -----------------------------------------------------------------------
% A trick: how to obtain the effect of the function 1(t) without applying
% the heaviside() MATLAB function
% (a) let's generate a logical array, containing "0" when the time instants
%     in tSET assume a value less than tDELAY, and "1" when the values in
%     tSET are greater than or equal to tDELAY
delayedSEQ = (tSET >=10);
% (b) exploiting the "automatic" casting feature 
%     (from "logical" to "double"), let's compute
%     the array containing the values of the 2nd input
u2SEQ = U2*delayedSEQ;
% -----------------------------------------------------------------------
 
% ------ assembling the input array -------
U  = [u1SEQ; u2SEQ];
% ------ assembling the input array -------
% ----------------- generation of the input sequences --------------------

Let's analyse the size and class of the variables we have used so far.

whos tSET u1SEQ delayedSEQ u2SEQ U
  Name            Size                Bytes  Class      Attributes

  U               2x20000            320000  double               
  delayedSEQ      1x20000             20000  logical              
  tSET            1x20000            160000  double               
  u1SEQ           1x20000            160000  double               
  u2SEQ           1x20000            160000  double               

Have we generated input sequences that have the correct behaviour? Let's examine the plots of the input signals.

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

plot(tSET, U, 'Linewidth', 2);

ylim([-3.5, +2.5]);

xlim([0,30]);

grid on;

ylabel('Inputs')

xlabel('Time $t$ [sec]', 'Interpreter','latex');

legend('$u_1(t)$','$u_2(t)$','Interpreter','latex', 'Location', 'best');

The LTI System

matA = [         0,            1,           0,           0;
        -(K1+K)/M1,  - (B1+B)/M1,        K/M1,        B/M1;
                 0,            0,           0,           1;
              K/M2,         B/M2,  -(K2+K)/M2,  -(B2+B)/M2];
 
matB = [    0,      0;
         1/M1,      0;
            0,      0;
            0,  -1/M2];
 
matC = [ 1,  0,  0,  0;
         0,  0,  1,  0];
 
matD = zeros(2, 2); % 2 outputs, 2 inputs
 
mechSYS = ss(matA, matB, matC, matD)
mechSYS =
 
  A = 
             x1        x2        x3        x4
   x1         0         1         0         0
   x2    -0.191   -0.0536     0.091    0.0036
   x3         0         0         0         1
   x4    0.0364   0.00144   -0.0764  -0.02144
 
  B = 
         u1    u2
   x1     0     0
   x2     1     0
   x3     0     0
   x4     0  -0.4
 
  C = 
       x1  x2  x3  x4
   y1   1   0   0   0
   y2   0   0   1   0
 
  D = 
       u1  u2
   y1   0   0
   y2   0   0
 
Continuous-time state-space model.
Model Properties

Free State and Output Movements

[YmechFREE, ~, XmechFREE] = initial(mechSYS, x0IC, tSET);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET, XmechFREE, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('State $x(t)$', 'Interpreter','latex');

title('Free State Movement')

legend('$x_1(t)$', '$x_2(t)$', '$x_3(t)$', '$x_4(t)$', 'Interpreter', 'latex', 'Location', 'best')

subplot(2, 1, 2)

plot(tSET, YmechFREE, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('Output $y(t)$', 'Interpreter','latex');

title('Free Output Movement')

legend('$y_1(t)$', '$y_2(t)$', 'Interpreter', 'latex', 'Location', 'best')

Forced State and Output Movements

[YmechLSIM, ~,XmechLSIM] = lsim(mechSYS, U, tSET);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET, XmechLSIM, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('State $x(t)$', 'Interpreter','latex');

title('Forced State Movement')

legend('$x_1(t)$', '$x_2(t)$', '$x_3(t)$', '$x_4(t)$', 'Interpreter', 'latex', 'Location', 'best')

subplot(2, 1, 2)

plot(tSET, YmechLSIM, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('Output $y(t)$', 'Interpreter','latex');

title('Forced Output Movement');

legend('$y_1(t)$', '$y_2(t)$', 'Interpreter', 'latex', 'Location', 'best')

Important Remark

Note the discontinuity in the first derivative of

at the time instant

sec. At that instant the second input starts influencing the behaviour of the system. You may note a discontinuity in the acceleration of the mass

(i.e., the first derivative of the speed of

State and Output Movements

To determine the overall movement (encompassing both free and forced movements) of the state variables and the overall movement of the outputs, we must sum the contributions from both free and forced movements.

XmechTOT = XmechFREE + XmechLSIM;

YmechTOT = YmechFREE + YmechLSIM;

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET, XmechTOT, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('State $x(t)$', 'Interpreter','latex');

title('State Movement')

legend('$x_1(t)$', '$x_2(t)$', '$x_3(t)$', '$x_4(t)$', 'Interpreter', 'latex', 'Location', 'best')

subplot(2, 1, 2)

plot(tSET, YmechTOT, 'Linewidth', 2);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('Output $y(t)$', 'Interpreter','latex');

title('Output Movement');

legend('$y_1(t)$', '$y_2(t)$', 'Interpreter', 'latex', 'Location', 'best')

Important Remark

Note the discontinuity in the first derivative of

at the time instant

sec. At that instant the second input starts influencing the behaviour of the system. You may note a discontinuity in the acceleration of the mass

(i.e., the first derivative of the speed of

Alternative Approach to Obtain Directly the Overall State and Output Movements

You may exploit the alternative syntax of the MATLAB command lsim to compute the overall state and output movements in one shot:

[YmechALT_TOT, ~,XmechALT_TOT] = lsim(mechSYS, U, tSET, x0IC);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET, XmechALT_TOT, 'Linestyle', '-', 'Linewidth', 2.5);

hold on;

plot(tSET, XmechTOT, 'Linestyle', '-.','Linewidth', 1.5);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('State $x(t)$', 'Interpreter','latex');

title('State Movement')

legend(' ALT $x_1(t)$', 'ALT $x_2(t)$', 'ALT $x_3(t)$', 'ALT $x_4(t)$',...

'$x_1(t)$', '$x_2(t)$', '$x_3(t)$', '$x_4(t)$','Interpreter', 'latex', 'Location', 'best')

subplot(2, 1, 2)

plot(tSET, YmechALT_TOT, 'Linestyle', '-', 'Linewidth', 2.5);

hold on;

plot(tSET, YmechTOT, 'Linestyle', '-.','Linewidth', 1.5);

grid on;

xlabel('Time $t$ [sec]', 'Interpreter','latex');

ylabel('Output $y(t)$', 'Interpreter','latex');

title('Output Movement');

legend('ALT $y_1(t)$', 'ALT $y_2(t)$', ...

'$y_1(t)$', '$y_2(t)$','Interpreter', 'latex', 'Location', 'best')

As you note, we obtained the same result by summing up free and forced movements or computing the overall movements directly.

A Linear Discrete Time System: a Supply Chain System

The dynamics of the system is given by the state-space model derived from conservation balances (of mass, energy, etc.).

Spurchases the quantity of raw material at each month k;
A fraction of raw material is discarded; a fraction is shipped to producer P;
A fraction of product is sold by P to retailer R, a fraction is discarded;
Retailer Rreturns a fraction of defective products every month, and sells a fraction to customers.

Paying attention to the material flow in each work step, we obtain the following description:

where

month counter
raw material in S
products in P
products in R
products sold to customers

The following constraints apply to the model parameters

and

Proposal: create a MATLAB model and a live script to analyse the evolution of this discrete time dynamical system. Use the following parameters to configure the system model:

choose random values for the system parameters , ensuring that constraints on the parameters are respected;
assuming that , and can take on non-integer values, choose for the initial state three random values in the range ;
input of raw material:
simulation time interval [months]:

A Solution Proposal

System Parameters

clear variables
 
% ------- the random number generator ----------------------------------
% Doing so, we are sure to obtain the same (random) sequence of values at
% every run of the live-script.
rng(481516);
% ------- the random number generator ----------------------------------
 
% ----- the system parameters alpha1, alpha2, and beta3 --------------------
alpha1 = rand; % a random number, uniformly distributed in the interval (0 1)
alpha2 = rand; % a random number, uniformly distributed in the interval (0 1)
beta3  = rand; % a random number, uniformly distributed in the interval (0 1)
 
% --- Now the remaining parameters ---
 
% --- a watchdog ----
MAXhits = 1000; 
% after MAXhits attempts, the watchdog breaks any loop,
% and it throws a warning message
 

% --- delta1 ---
delta1 = rand; 
hit = 0; % initialise the watchdog counter
 
while (alpha1+delta1 > 1)
    delta1 = rand;
    hit = hit+1;
    if (hit > MAXhits)
        error('Delta1 loop failed! Reached the MAX number of attempts! Aborting! ')
    end
end
% leaving the "while" loop we are sure that the constraint
% 0 < alpha1+delta1 <=1
% has been met
ad1 = alpha1 + delta1;
assert( (ad1>0)&&(ad1<=1) )
 

% --- delta2 ---
delta2 = rand; 
hit = 0; % initialise the watchdog counter
 
while (alpha2+delta2 > 1)
    delta2 = rand;
    hit = hit+1;
    if (hit > MAXhits)
        error('Delta2 loop failed! Reached the MAX number of attempts! Aborting! ')
    end
end
% leaving the "while" loop we are sure that the constraint
% 0 < alpha2+delta2 <=1
% has been met
ad2 = alpha2 + delta2;
assert( (ad2>0)&&(ad2<=1) )

% --- gamma3 ---
gamma3 = rand; 
hit = 0; % initialise the watchdog counter
 
while (beta3+gamma3 > 1)
    gamma3 = rand;
    hit = hit+1;
    if (hit > MAXhits)
        error('Gamma3 loop failed! Reached the MAX number of attempts! Aborting! ')
    end
end
% leaving the "while" loop we are sure that the constraint
% 0 < alpha2+gamma3 <=1
% has been met
ad3 = beta3 + gamma3;
assert( (ad3>0)&&(ad3<=1) )

% ---------- the initial state --------------
% x1 <--> raw material in S
% x2 <--> products in P
% x3 <--> products in R
min_X_VAL = 100;
max_X_VAL = 500;
x0IC_supplyChain = (max_X_VAL - min_X_VAL) * rand([3,1]) + min_X_VAL .* ones([3,1]);
% ---------- the initial state --------------

Important Remark: Linear Transformations and Random Number Generators

The code on line 349 implements a linear transformation to alter the range of the uniformly distributed random generator rand().

In fact, rand() provides pseudorandom numbers in the interval . Refer to the documentation for details.
To obtain uniformly distributed random numbers in the general interval , we must map the uniformly distributed random numbers (values between 0 and 1) to the different range . The most straightforward approach consists of applying a linear transformation, as follows

% ----------- simulation time interval -------------------------------
% Let's create an array of time instants. We will determine 
% both the free state movement (the zero-input movement) and the forced
% state movement (the zero-state movement) in these time instants.
 
Tstart = 0; % the starting month of the simulation
Tstop = 36; % the simulation must end when the Tstop months expire.
MonthlyStep = 1; % the time period corresponds to one month
tSET_supplyChain = (Tstart:MonthlyStep:Tstop); 
% The array containing the time instants
% at which the values of free and forced
% movements will be determined.
% ----------- simulation time interval -------------------------------

Input Generation

Let us assume that

can take non-integer values. Moreover, we can avoid to use the Heaviside function in the computation

uSTEPterm = 75.0;

uSINterm = 35.0;

uANG_FREQ = 6;

uSEQsupplyChain = uSTEPterm * ones(size(tSET_supplyChain)) + uSINterm*sin(uANG_FREQ * tSET_supplyChain);

whos uSEQ

figure('Units', 'pixels', 'Position',[1, 1, 1280, 980]);

plot(tSET_supplyChain, uSEQsupplyChain, 'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 30);

grid on;

xlabel('Month'); ylabel('Raw Material')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)]);

ylim([0.8*min(uSEQsupplyChain), 1.05*max(uSEQsupplyChain)]);

The LTI System

matA = [ (1-alpha1-delta1),                  0,                 0;
                    alpha1,  (1-alpha2-delta2),             beta3;
                         0,             alpha2,   (1-beta3-gamma3)];
 
matB = [ 1;
         0;
         0 ];
 
matC = [ 0,    0,  gamma3];
 
matD = 0;
 
% MonthlyStep = 1; % the time period corresponds to one month
% NOTE" refer to the ss documentation:
% https://it.mathworks.com/help/releases/R2025a/control/ref/ss.html?lang=en#d126e261479
SupplyChainSYS = ss(matA, matB, matC, matD, MonthlyStep, 'TimeUnit', 'months')
SupplyChainSYS =
 
  A = 
            x1       x2       x3
   x1  0.04515        0        0
   x2  0.05612   0.2093   0.2848
   x3        0   0.6405   0.4806
 
  B = 
       u1
   x1   1
   x2   0
   x3   0
 
  C = 
           x1      x2      x3
   y1       0       0  0.2346
 
  D = 
       u1
   y1   0
 
Sample time: 1 months
Discrete-time state-space model.
Model Properties

Free State and Output Movements

[YsupplyFREE, ~, XsupplyFREE] = initial(SupplyChainSYS, x0IC_supplyChain, tSET_supplyChain);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_supplyChain, XsupplyFREE, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('Free State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

subplot(2, 1, 2)

plot(tSET_supplyChain, YsupplyFREE, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Free Output Movement')

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

Forced State and Output Movements

[YsupplyLSIM, ~,XsupplyLSIM] = lsim(SupplyChainSYS, uSEQsupplyChain, tSET_supplyChain);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_supplyChain, XsupplyLSIM, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('Forced State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

subplot(2, 1, 2)

plot(tSET_supplyChain, YsupplyLSIM, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Forced Output Movement');

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

State and Output Movements

XsupplyTOT = XsupplyFREE + XsupplyLSIM;

YsupplyTOT = YsupplyFREE + YsupplyLSIM;

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_supplyChain, XsupplyTOT, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

subplot(2, 1, 2)

plot(tSET_supplyChain, YsupplyTOT, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time $k$ [month]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Output Movement');

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_supplyChain(1), tSET_supplyChain(end)])

A Linear Discrete Time System: a linear Finite Difference Equation

The state-space model is derived from a difference equation model (an Input/Output model):

Letting

one gets:

Proposal: create a MATLAB model and a live script to analyse the evolution of this discrete time dynamical system. Use the following parameters to configure the system model:

choose random values for the initial state
assuming as input the sequence:
simulation time interval:

A Solution Proposal

System Parameters

clear variables
 
% ------- the random number generator ----------------------------------
% Doing so, we are sure to obtain the same (random) sequence of values at
% every run of the live-script.
rng(162342);
% ------- the random number generator ----------------------------------
 
% ---------- the initial state --------------
% x1(k)  <--> y(k-3)
% x2(k)  <--> y(k-2)
% x3(k)  <--> y(k-1)
 
x0IC_diffEQ = rand([3,1]);
% ---------- the initial state --------------

% ----------- simulation time interval -------------------------------
% Let's create an array of time instants. We will determine 
% both the free state movement (the zero-input movement) and the forced
% state movement (the zero-state movement) in these time instants.
 
Tstart = 0; % the starting month of the simulation
Tstop = 25; % the simulation must end when the Tstop months expire.
timeStep = 1; % the time period corresponds to one step forward for the solution of the difference equation
tSET_diffEQ = (Tstart:timeStep:Tstop); 
% The array containing the time instants
% at which the values of free and forced
% movements will be determined.
% ----------- simulation time interval -------------------------------

Input Generation

Let us assume that

can take non-integer values. Moreover, we can avoid to use the Heaviside function in the computation

uSTEPterm = 1.25;

uPOWterm = -2/3;

uSEQdiffEQ = uSTEPterm * ones(size(tSET_diffEQ)) + uPOWterm .^ (0:numel(tSET_diffEQ)-1);

whos uSEQ

figure('Units', 'pixels', 'Position',[1, 1, 1280, 980]);

plot(tSET_diffEQ, uSEQdiffEQ, 'LineStyle', 'none', 'Marker', '.', 'MarkerSize', 30);

grid on;

xlabel('Time Instant k'); ylabel('Input u(k)')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)]);

ylim([0.8*min(uSEQdiffEQ), 1.05*max(uSEQdiffEQ)]);

The LTI System

matA = [ 0,   1,  0;
         0,   0,  1;
         1,  -2,  3 ];
 
matB = [ 0;
         0;
         6 ];
 
matC = [ 1,  -2,  3];
 
matD = 6;
 
% timeStep = 1; % the time period corresponds to one step forward for the solution of the difference equation
unspecifiedTs = -1; % it's a discrete time system with un unspecified sampling time
diffEQ_Sys = ss(matA, matB, matC, matD, unspecifiedTs )
diffEQ_Sys =
 
  A = 
       x1  x2  x3
   x1   0   1   0
   x2   0   0   1
   x3   1  -2   3
 
  B = 
       u1
   x1   0
   x2   0
   x3   6
 
  C = 
       x1  x2  x3
   y1   1  -2   3
 
  D = 
       u1
   y1   6
 
Sample time: unspecified
Discrete-time state-space model.
Model Properties

Free State and Output Movements

[YdiffFREE, ~, XdiffFREE] = initial(diffEQ_Sys, x0IC_diffEQ, tSET_diffEQ);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_diffEQ, XdiffFREE, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('Free State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

subplot(2, 1, 2)

plot(tSET_diffEQ, YdiffFREE, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Free Output Movement')

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

Forced State and Output Movements

[YdiffLSIM, ~,XdiffLSIM] = lsim(diffEQ_Sys, uSEQdiffEQ, tSET_diffEQ);

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_diffEQ, XdiffLSIM, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('Forced State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

subplot(2, 1, 2)

plot(tSET_diffEQ, YdiffLSIM, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Forced Output Movement');

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

State and Output Movements

XdiffTOT = XdiffFREE + XdiffLSIM;

YdiffTOT = YdiffFREE + YdiffLSIM;

figure('Units', 'pixels', 'Position', [1,1,1920, 1080]);

subplot(2, 1, 1)

plot(tSET_diffEQ, XdiffTOT, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('State $x(k)$', 'Interpreter','latex');

title('State Movement')

legend('$x_1(k)$', '$x_2(k)$', '$x_3(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

subplot(2, 1, 2)

plot(tSET_diffEQ, YdiffTOT, 'LineStyle', 'none', 'Marker','.', 'MarkerSize', 20);

grid on;

xlabel('Time instant $k$ [ ]', 'Interpreter','latex');

ylabel('Output $y(k)$', 'Interpreter','latex');

title('Output Movement');

legend('$y_1(k)$', 'Interpreter', 'latex', 'Location', 'best')

xlim([tSET_diffEQ(1), tSET_diffEQ(end)])

A Nonlinear Continuous-Time Dynamic System: a Double Pendulum

The double pendulum consists of two pendulums and two rigid rods. The pendulums are attached to each other, as in the following figure

Assumptions regarding the dynamic system's model:

the pendulum rods are mass-less and rigid;
the pendulum masses are point masses;
the only force acting on pendulums is the gravity force.

The Model

The ODEs describing the system's motion can be derived using the direct Newtonian method (refer to [3], for example):

where

and are the angles of the first and latter pendulum, respectively;
and are the (constant) lengths of the two rods;
and are the (constant) masses of the pendulums;
g is the gravitational constant.

Note: for simplicity, it is assumed that there is no input. Gravity acceleration is NOT considered an input of the system but simply a parameter of the ODEs.

The dynamic model of the double pendulum is a set of two second-order nonlinear Ordinary Differential Equations (ODEs).

As outputs of the system consider simply:

System's Parameters

clear variables
 
g = 9.80665;  % Acceleration due to gravity (m/s^2)
 
% first set of parameters
% L1 = 1.0;  % Length of the first pendulum arm (m)
% L2 = 1.0;  % Length of the second pendulum arm (m)
% m1 = 1.0;  % Mass of the first pendulum bob [= the ball, attached to the end of the rod]  (kg)
% m2 = 1.0;  % Mass of the second pendulum bob (kg)
 
% a different set
L1 = 5.0;  % Length of the first pendulum arm (m)
L2 = 2.5;  % Length of the second pendulum arm (m)
m1 = 1.0;  % Mass of the first pendulum bob [= the ball, attached to the end of the rod]  (kg)
m2 = 2.0;  % Mass of the second pendulum bob (kg)

How to Solve ODEs in MATLAB

Critical Remark: from an N-th Order ODEs Set to a Set of First-Order ODEs

The order of the ODE equations is crucial: to be able to solve numerically in MATLAB a system of (linear or nonlinear) ODEs of order n (with

), we must always rewrite the n-th order ODEs set as a set of first-order ODEs, properly coupled each to the others. Once we have obtained the set of first-order ODEs, we can

describe this set using the MATLAB command ode, in case the ODEs are nonlinear and then solve the ODEs;
exploit the Control System Toolbox functions (as, for example, lsim) in case the ODEs are describing a linear time-invariant system (refer to the LTI Systems devoted live script in this collection).

The MATLAB Command ode

For help on using the command ode, and examples, refer to the online documentation or execute the following piece of code in the MATLAB Command Window

doc ode

Briefly, the ode command allows us to solve Initial Value Problems (IVPs) of the form

where y is a (vector) function (the solution of the ODE we are looking for) and p is a vector of parameters.

Thus, since the ODEs

are already a set of first-order ODEs (nonlinear ODEs), we can describe the ODEs set structure using the command ode.

Important Remark on Notation

It's not mandatory to obey the notation

. Indeed, we will use a different notation for IPVs, when they describe the evolution of a dynamic system (linear or nonlinear), keeping the notation x for the state variables (the variables affected by the ODEs), and using the notation y for the system's outputs. Thus, our default notation for an IVP is

How to Obtain a Set of First-Order ODEs

Given the ODEs set

, let's transform it into a set of first-order ODEs with proper initial conditions:

The dynamic system is mechanical; thus, the positions and velocities (linear, angular) of rigid bodies’ centers of mass are candidates for state variables.
The rods are mass-less, so the state variables are the angular positions and velocities of the two masses.

Using this set of variables, let's transform the Eq.s

into a set of first-order ODEs, i.e. a state-space description

The IVP, we would like to solve in MATLAB, is described by the Eqs

together with the initial conditions

where

is the initial time instant.

The IVP Initial Conditions

t0    =  0; % starting time instant [s]
tStop = 20; % end of the simulation [s]
 
%
barTheta1 = pi/2; % Initial angle of the first pendulum (radians)
barOmega1 = 0;    % Initial angular velocity of the first pendulum (rad/s)
 
barTheta2 = pi/2; % Initial angle of the second pendulum (radians)
barOmega2 = 0;    % Initial angular velocity of the second pendulum (rad/s)
%}
 
%{
% a different IC configuration
barTheta1 = pi/2; % Initial angle of the first pendulum (radians)
barOmega1 = 0;    % Initial angular velocity of the first pendulum (rad/s)
 
barTheta2 = deg2rad(-10); % Initial angle of the second pendulum (radians) *** try -8 deg ***
barOmega2 = 0;    % Initial angular velocity of the second pendulum (rad/s)
%}

A Solution Proposal

Describing the ODEs in MATLAB

%
DoublePend_ODE = ode;                               % let's initialise an ODE object
DoublePend_ODE.Parameters   = [g, L1, L2, m1, m2];  % the parameters of the differential equation
% a set of parameters
% L1  Length of the first pendulum arm (m)
% L2  Length of the second pendulum arm (m)
% m1  Mass of the first pendulum bob [= the ball, attached to the end of the rod]  (kg)
% m2  Mass of the second pendulum bob (kg)
 
DoublePend_ODE.InitialValue = [barTheta1, barOmega1, barTheta2, barOmega2]; % the initial conditions
DoublePend_ODE.InitialTime = t0;  % starting time instant
                                  % the initial time instant of the solution, where 
                                  % the initial conditions are applied
DoublePend_ODE.ODEFcn = @doublePEND_ODE; % doublePEND_ODE is the MATLAB function describing the ODE, 
                                         % we would like to solve. 
                                         % The code of doublePEND_ODE is at the end
                                         % of this live script.
 

Simulation of the System - Solving the ODE

Given the MATLAB description of the IVP we would solve, let's determine a numeric solution by applying the MATLAB command solve.

For help on using the command solve, and examples, refer to the online documentation or execute the following piece of code in the MATLAB Command Window

doc solve

Briefly, the solve command computes the solution of the considered IVP F in the time interval [t0, tf], returning the solution S at each time step chosen by the solver:

S = solve(F, t0, tf);

S  = solve(DoublePend_ODE, t0, tStop); % solve the IPV described by F_ODE 
                           % in the time interval [t0, tf]
 
doubleP_x = S.Solution; % the result of the IPv: the state movement
t_stamps = S.Time;   % time instants corresponding to the solution samples 
 
doubleP_x1 = doubleP_x(1,:);
doubleP_x2 = doubleP_x(2,:);
doubleP_x3 = doubleP_x(3,:);
doubleP_x4 = doubleP_x(4,:);
 
y1 = doubleP_x1; % the model output y_1
y2 = doubleP_x3; % the model output y_2
 
tMIN = min(t_stamps); tMAX = max(t_stamps);

Retrieving and Showing the Result

figure('Name','Double Pendulum Simulation','Units','normalized','Position',[0.1, 0.1, 0.75, 0.95]);

h1 = subplot(2,1,1);

plot(t_stamps, doubleP_x1, 'LineWidth',2); % first state variable

hold on

plot(t_stamps, doubleP_x2, 'LineWidth',2); % second state variable

plot(t_stamps, doubleP_x3, 'LineWidth',2); % third state variable

plot(t_stamps, doubleP_x4, 'LineWidth',2); % fourth state variable

xlabel('time [h]'); ylabel('state variables'); % labels on the plot axes

grid on; currAX = axis(gca); currAX(1) = tMIN; currAX(2) = tMAX;

axis(currAX);set(gca, 'FontSize', 14);

legend('x_1(t)', 'x_2(t)', 'x_3(t)', 'x_4(t)', 'Location', 'best')

h2 = subplot(2,1,2);

plot(t_stamps, y1, 'LineWidth',2); % double pendulum outputs y_1 vs time

hold on

plot(t_stamps, y2, 'LineWidth',2); % double pendulum outputs y_2 vs time

xlabel('time [h]'); ylabel('ang. position of the pendulum bobs [rad]'); % labels on the plot axes

grid on;

currAX = axis(gca); currAX(1) = tMIN; currAX(2) = tMAX;

axis(currAX);set(gca, 'FontSize', 14);

legend('y_1(t)', 'y_2(t)', 'Location', 'best')

linkaxes([h1 h2], 'x');

The Double Pendulum Model: The ODE as MATLAB Function

Let's define as a MATLAB function the first-order ODE system corresponding to

The MATLAB function must provide the first derivative vector

, given the values of the variables

and

at the current time instant t and the parameters' values.

% Define the function for the double pendulum ODEs
function dxdt = doublePEND_ODE(t, x, p)
 
    g  = p(1);
    L1 = p(2);
    L2 = p(3);
    m1 = p(4);
    m2 = p(5);
    % Unpack the state variables
    theta1 = x(1);
    omega1 = x(2);
    theta2 = x(3);
    omega2 = x(4);
 
    % Equations of motion for the double pendulum
    dxdt = zeros(4, 1);
    dxdt(1) = omega1;
    dxdt(2) = (-g * (2 * m1 + m2) * sin(theta1) - m2 * g * sin(theta1 - 2 * theta2) - 2 * sin(theta1 - theta2) * m2 * (omega2^2 * L2 + omega1^2 * L1 * cos(theta1 - theta2))) / (L1 * (2 * m1 + m2 - m2 * cos(2 * theta1 - 2 * theta2)));
    dxdt(3) = omega2;
    dxdt(4) = (2 * sin(theta1 - theta2) * (omega1^2 * L1 * (m1 + m2) + g * (m1 + m2) * cos(theta1) + omega2^2 * L2 * m2 * cos(theta1 - theta2))) / (L2 * (2 * m1 + m2 - m2 * cos(2 * theta1 - 2 * theta2)));
end
 