State Space Descriptions: Criteria and Examples

Table of Contents

Introduction Modeling Based on Conservation Balances A Nonlinear Continuous Time System: a Continuous Bio-Reactor Cleaning the MATLAB Workspace Assigning the System Parameters Initial Condition and Input Configuration Simulation Time Interval Solving an Ordinary Differential Equation in MATLAB Important Remark The MATLAB Command ode Important Remark on Notation Simulation of the System - Solving the ODE Retrieving and Showing the Result Modeling Based on Input-Output Mathematical Description A Nonlinear Continuous Time System: a Second Order Ordinary Differential Equation Parameters Initial values Input Cleaning the MATLAB Workspace Assigning the System Parameters Initial Condition and Input Configuration Simulation Time Interval Solving the ODE in MATLAB The MATLAB Command ode Important Remark on Notation Simulation of the System - Solving the ODE Retrieving and Showing the Result A Nonlinear Discrete Time System: a Third Order Difference Equation Cleaning the Matlab Workspace Assigning the Initial Condition and Input Configuration Simulation Time Interval The Input Sequence The Model and the Simulation of the System Retrieving and Showing the Result A Proposal The Bioreactor Model: The ODE as MATLAB Function The Nonlinear Input-Output Model: The ODE as MATLAB Function The Third Order Nonlinear Difference Equation: the Difference Equation as a MATLAB Function

Introduction

This live script uses MATLAB to implement a state space description for some dynamic (linear or nonlinear, time-invariant or time-variant, continuous time) systems.

What you will learn:

Given a state space description of a dynamic system, how to build a MATLAB model for such a system;
Tune a MATLAB model, setting the model parameters (or the simulation parameters) using the command ode;
How to simulate a MATLAB model of a dynamic system, and retrieve the results using the MATLAB command solve.

Modeling Based on Conservation Balances

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

A Nonlinear Continuous Time System: a Continuous Bio-Reactor

Example of fermenter-bioreactor (Source: LAMBDA Instruments)

The dynamics of the process are derived from the conservation balances for the biomass with concentration

and substrate mass with concentration

where

is the biomass concentration [g/l];
is the substrate concentration [g/l];
u is the feed flow rate [l/h]

and the parameters assume the following values

volume V:
substrate feed concentration :
yield coefficient Y:
maximal grow rate :
saturation parameter :
inhibition parameter :

Finally, the initial condition for the bioreactor corresponds to the following values for the state variables:

and the input

is defined as follows

where

is the Heaviside function

Cleaning the MATLAB Workspace

close all
clear
clc
 

Assigning the System Parameters

V = 4;     % volume [l]
S_F = 10;  % substrate feed concentration [g/l]
Y = 0.5;   % yield coefficient
muMAX = 1; % maximal grow rate [l/h]
K1 = 0.03; % saturation parameter [g/l]
K2 = 0.5;  % inhibition parameter [g/l]

Initial Condition and Input Configuration

x2_0 = -(K1 - (K1*(K2*S_F^2 + S_F + K1))^(1/2))/(K2*S_F + 1); % initial substrate concentration [g/l]
x1_0 =(S_F-x2_0)*Y;                                           % initial biomass concentration [g/l]
 
U_M = V*muMAX*x2_0/(K2*x2_0^2+x2_0+K1);   % mean input flow rate [l/h]
                                     
u_A = 0.75e-3;    % amplitude of the sinusoidal wave input component [l/h]
T_u = 2.5 * 3600; % time period of the periodic input component [seconds]
 

Simulation Time Interval

T_start = 0.0;        % starting time instant
T_stop = 10.0 * 3600; % end of the simulation after 10 h

Solving an Ordinary Differential Equation in MATLAB

The dynamic model of the bioreactor is a system of two first-order Ordinary Differential Equations (ODEs).

where

is the Heaviside function

Important Remark

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

), we must always rewrite the n-th order ODEs set as a system 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.

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

Putting the ODEs in MATLAB: Let's define the system of first-order ODEs with given initial conditions in MATLAB, using the command ode

where

is the Heaviside function

%
BioR_ODE = ode;                                  % let's initialise an ODE object
BioR_ODE.Parameters   = [V, S_F, Y, muMAX, ...
                         K1, K2, U_M, u_A, T_u];  % the parameters of the differential equation
% V = 4;     % volume [l]
% S_F = 10;  % substrate feed concentration [g/l]
% Y = 0.5;   % yield coefficient
% muMAX = 1; % maximal grow rate [l/h]
% K1 = 0.03; % saturation parameter [g/l]
% K2 = 0.5;  % inhibition parameter [g/l]
 
BioR_ODE.InitialValue = [x1_0, x2_0];  % the initial conditions
BioR_ODE.InitialTime = T_start;        % starting time instant
                                       % the initial time instant of the solution, where 
                                       % the initial conditions are applied
BioR_ODE.ODEFcn = @BioreactorODE;  % BioreactorODE is the MATLAB function describing the ODE, 
                                % we would like to solve. 
                                % The code of BioreactorODE 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(BioR_ODE, T_start, T_stop); % solve the IVP described by F_ODE 
                           % in the time interval [t0, tf]
 
bioreactor_x = S.Solution; % the result of the IVP
t_stamps = S.Time;   % time instants corresponding to the solution samples 
 
bioreactor_x1 = bioreactor_x(1,:);
bioreactor_x2 = bioreactor_x(2,:);
 
y1 = bioreactor_x1; % the model output y_1
y2 = bioreactor_x2; % the model output y_2
u = U_M+u_A*sin((2*pi)*(t_stamps/T_u));  % the input u <-- NOTE: We are forced to recalculate the input 
 
t_stamps = t_stamps/3600; % time instants [h]
% time instants corresponding to the solution samples 

Retrieving and Showing the Result

figure('Name','Bioreactor Simulation','Units','pixels','Position',[1, 1, 1280, 1080]);

h1 = subplot(3,1,1);

plot(t_stamps, y1, 'LineWidth',2, 'Color',[0.9290 0.6940 0.1250]); % bioreactor outputs y_1 vs time

xlabel('time [h]'); ylabel('biomass concentration [g/l]'); % labels on the plot axes

grid on;

tMIN = min(t_stamps); tMAX = max(t_stamps);

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

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

h2 = subplot(3,1,2);

plot(t_stamps, y2, 'LineWidth',2, 'Color',[0.8500 0.3250 0.0980]); % bioreactor outputs y_2 vs time

xlabel('time [h]'); ylabel('substrate concentration [g/l]'); % labels on the plot axes

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

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

h3 = subplot(3,1,3);

plot(t_stamps, u, 'LineWidth',2, 'Color',[0.4660 0.6740 0.1880]); % bioreactor outputs y_2 vs time

xlabel('time [h]'); ylabel('feed flow rate (input) [g/l]'); % labels on the plot axes

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

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

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

Modeling Based on Input-Output Mathematical Description

The state-space model is derived from an input/output model based on differential equations or from an in-out difference equation model.

A Nonlinear Continuous Time System: a Second Order Ordinary Differential Equation

Let us consider the nonlinear dynamical system defined by the differential equation

This model can be expressed in state-space form as

where

by assuming

Parameters

Assume the following values for the differential equation parameters:

Initial values

Moreover, the initial values for the state variables are

Input

Apply as input

the following signal:

where

is the Heaviside function

Cleaning the MATLAB Workspace

close all
clear
clc
 

Assigning the System Parameters

alpha1 = 3.0;  %#ok<*NASGU> 
alpha0 = 2.0;  
gamma1 = 0.5;  
beta0 = 2.0; 
gamma2 = 0.25; 

Initial Condition and Input Configuration

x1_0 =  +2.50; % initial value of the 1st state variable
x2_0 = -24.48; % initial value of the 2nd state variable
 
u_expT = 1/4;  % time constant of the exponential input component
u_sinA = 2.50; % amplitude of the sinusoidal input component
u_sinO = 3.0;  % angular frequency of the input sinusoidal component [rad/s]
u_cosA = 7.25; % amplitude of the cosinusoidal input component
u_cosO = 2.0;  % angular frequency of the input cosinusoidal component [rad/s]
 

Simulation Time Interval

T_start2 = 0.0; % starting time instant [s]
T_stop2 = 5.0; % end of the simulation  [s]

Solving the ODE in MATLAB

By exploiting the definitions

we already transformed the second-order ODE

into the system of two first-order Ordinary Differential Equations (ODEs).

where

is the Heaviside function

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

Putting the ODEs in MATLAB: Let's define the system of first-order ODEs with given initial conditions in MATLAB, using the command ode

where

is the Heaviside function

%
IO_ODE = ode;                                  % let's initialise an ODE object
IO_ODE.Parameters   = [alpha1, alpha0, gamma1, beta0, ...
                         gamma2, u_expT, u_sinA, u_sinO, ...
                         u_cosA, u_cosO];  % the parameters of the differential equation
% alpha1 = 3.0;  
% alpha0 = 2.0;  
% gamma1 = 0.5;  
% beta0 = 2.0; 
% gamma2 = 0.25; 
% 
% u_expT = 1/4;  % time constant of the exponential input component
% u_sinA = 2.50; % amplitude of the sinusoidal input component
% u_sinO = 3.0;  % angular frequency of the input sinusoidal component [rad/s]
% u_cosA = 7.25; % amplitude of the cosinusoidal input component
% u_cosO = 2.0;  % angular frequency of the input cosinusoidal component [rad/s]
 
 
 
IO_ODE.InitialValue = [x1_0, x2_0];  % the initial conditions
% x1_0 =  +2.50; % initial value of the 1st state variable
% x2_0 = -24.48; % initial value of the 2nd state variable
 
IO_ODE.InitialTime = T_start2;        % starting time instant
                                     % the initial time instant of the solution, where 
                                     % the initial conditions are applied
% T_start = 0.0; % starting time instant [s]
% T_stop = 25.0; % end of the simulation  [s]
 
IO_ODE.ODEFcn = @InOutODE;  % InOutODE is the MATLAB function describing the ODE, 
                            % we would like to solve. 
                            % The code of InOutODE 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(IO_ODE, T_start2, T_stop2); % solve the IVP described by F_ODE 
                           % in the time interval [t0, tf]
 
IO_x = S.Solution; % the result of the IVP
t_stamps = S.Time;   % time instants corresponding to the solution samples 
 
IO_x1 = IO_x(1,:);
IO_x2 = IO_x(2,:);
 
y1 = IO_x1; % the model output y_1
 
u = exp(-t_stamps/u_expT).*(u_sinA.*sin(u_sinO.*t_stamps)+...
    u_cosA.*cos(u_cosO.*t_stamps));  % the input u <-- NOTE: We are forced to recalculate the input 
 

Retrieving and Showing the Result

figure('Name','State-Space from ODE','Units','pixels','Position',[1, 1, 1280, 1280]);

h1 = subplot(3,1,1);

plot(t_stamps, y1, 'LineWidth',2, 'Color',[0.9290 0.6940 0.1250]); % ODE output vs time

xlabel('time [s]'); ylabel('ODE output y(t)'); % labels on the plot axes

grid on;

tMIN = min(t_stamps); tMAX = max(t_stamps);

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

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

h2 = subplot(3,1,2);

plot(t_stamps, u, 'LineWidth',2, 'Color',[0.4660 0.6740 0.1880]); % bioreactor outputs y_2 vs time

xlabel('time [s]'); ylabel('ODE input u(t)'); % labels on the plot axes

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

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

h3 = subplot(3,1,3);

plot(t_stamps, IO_x1, 'LineWidth',2, 'Color',[0.8500 0.3250 0.0980]); % state variable x_1

hold on;

plot(t_stamps, IO_x2, 'LineWidth',2, 'Color',[0.3010 0.7450 0.9330]); % state variable x_2

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

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

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

legend('x_1(t)', 'x_2(t)', 'Location','best');

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

A Nonlinear Discrete Time System: a Third Order Difference Equation

Let us consider the nonlinear dynamical system defined by the difference equation

with initial conditions

and the input sequence

This model can be expressed in state-space form as

where

Cleaning the Matlab Workspace

close all
clear
clc
 

Assigning the Initial Condition and Input Configuration

x1_0 = +2.10; % initial value of the 1st state variable x1(k) <--> y(k-2)
x2_0 = -3.40; % initial value of the 2nd state variable x2(k) <--> y(k-1)
x3_0 = +0.25; % initial value of the 3rd state variable x3(k) <--> y(k)
 
u_Om1 = 3.00; % angular frequency of the 1st input sinusoidal component [rad/sample]
u_Om2 = 0.25; % angular frequency of the 2nd input sinusoidal component [rad/sample]
u_Om3 = 40.0; % angular frequency of the 3rd input sinusoidal component [rad/sample]
 

Simulation Time Interval

T_start = 0.0;   % starting time instant [s]
T_stop = 1250.0; % end of the simulation  [s]
Ts = 1;          % discrete time period [s]
timeT = (T_start:Ts:T_stop); % the finite sequence of simulation time instants

The Input Sequence

uK = sin(u_Om1 .* timeT) + sin(u_Om2 .* timeT) - cos(u_Om3 .*timeT); 
% timeT is starting from t=0, thus we DO NOT need to multiply uK by the term 
% heaviside(timeT) (this term is always equal to 1)

The Model and the Simulation of the System

This time, we create a MATLAB function that implements the state equations describing the original nonlinear difference equation and use it to compute the state and output discrete-time evolution.

x0 = [x1_0; x2_0; x3_0]; % the state IC
u_in = uK;               % the input sequence
 
[y_out, x_out] = Solve_ORD3_NONLIN_DifferenceEquation(timeT, x0, u_in );
 
x_1 = x_out(:,1); % 1st state variable: x_1(k) <--> y(k-2)
x_2 = x_out(:,2); % 2nd state variable: x_2(k) <--> y(k-1)
x_3 = x_out(:,3); % 3rd state variable: x_3(k) <--> y(k)

Retrieving and Showing the Result

figure('Name','State-Space from Difference Equations','Units','pixels','Position',[1, 1, 1280, 1280]);

h1 = subplot(3,1,1);

plot(timeT, y_out, 'lineStyle', ' -', 'LineWidth', 1,...

'Marker', 'o', 'Color',[0.9290 0.6940 0.1250], ...

'MarkerSize', 4, 'MarkerFaceColor',[0.9290 0.6940 0.1250],...

'MarkerEdgeColor', [0.9290 0.6940 0.1250]); % Diff. Eq output vs time

xlabel('time [s]'); ylabel('DE output y(k)'); % labels on the plot axes

grid on; zoom on

tMIN = min(timeT); tMAX = max(timeT);

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

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

h2 = subplot(3,1,2);

plot(timeT, u_in, 'LineStyle', '-.', 'LineWidth',1, 'Color',[0.4660 0.6740 0.1880], ...

'Marker', 'diamond', 'MarkerSize',4,'MarkerFaceColor',[0.4660 0.6740 0.1880],...

'MarkerEdgeColor',[0.4660 0.6740 0.1880]); % bioreactor outputs y_2 vs time

xlabel('time [s]'); ylabel('DE input u(k)'); % labels on the plot axes

grid on; zoom on

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

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

h3 = subplot(3,1,3);

plot(timeT, x_1, 'LineStyle', ':', 'LineWidth', 1,'Color',[0.8500 0.3250 0.0980], ...

'Marker','square','MarkerSize', 4, 'MarkerFaceColor', [0.8500 0.3250 0.0980],...

'MarkerEdgeColor',[0.8500 0.3250 0.0980]); % state variable x_1

hold on;

plot(timeT, x_2,'LineStyle', ':', 'LineWidth', 1,'Color',[0.3010 0.7450 0.9330],...

'Marker','pentagram','MarkerSize', 4, 'MarkerFaceColor', [0.3010 0.7450 0.9330],...

'MarkerEdgeColor', [0.3010 0.7450 0.9330]); % state variable x_2

plot(timeT, x_3,'LineStyle', ':', 'LineWidth', 1,'Color',[0.4940 0.1840 0.5560],...

'Marker', '^','MarkerSize', 4, 'MarkerFaceColor', [0.4940 0.1840 0.5560],...

'MarkerEdgeColor', [0.4940 0.1840 0.5560]); % state variable x_2

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

grid on; zoom on

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

axis(currAX);

legend('x_1(k)', 'x_2(k)', 'x_3(k)', 'Location','best');

linkaxes([h1 h2 h3], 'x');set(gca, 'FontSize', 10);

A Proposal

Using what is shown in this live script, create the Simulink models and write the Matlab code necessary to configure them, run the simulation and retrieve the results, for the examples proposed in the live script named HandsON_StateSpaceDescription.m

The Bioreactor 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.

 
function dxdt = BioreactorODE(t, x, p)
    % BioR_ODE.Parameters   = [V, S_F, Y, muMAX, ...
    %                      K1, K2, U_M, u_A, T_u];  % the parameters of the differential equation
    % retrieving the parameters
    V_val     = p(1);
    SF_val    = p(2);
    Y_val     = p(3);
    muMAX_val = p(4);
    K1_val    = p(5);
    K2_val    = p(6);
    UM_val    = p(7);
    uA_val    = p(8);
    TU_val    = p(9);
 
    muVAL = muMAX_val*x(2)/(K2_val*(x(2)^2) + x(2) + K1_val);
    ut    = UM_val+uA_val*sin((2*pi)*(t/TU_val));
    dxdt = [muVAL*x(1)-(x(1)*ut)/V_val; ...
            -(muVAL*x(1))/Y_val+((SF_val-x(2))*ut)/V_val];
 
end % function

The Nonlinear Input-Output 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.

function dxdt = InOutODE(t, x, p)
   % IO_ODE.Parameters   = [alpha1, alpha0, gamma1, beta0, ...
   %                       gamma2, u_expT, u_sinA, u_sinO, ...
   %                       u_cosA, u_cosO];  % the parameters of the differential equation
   %  % retrieving the parameters
    alpha1_val = p(1);
    alpha0_val = p(2);
    gamma1_val = p(3);
    beta0_val  = p(4);
    gamma2_val = p(5);
    u_expT_val = p(6);
    u_sinA_val = p(7);
    u_sinO_val = p(8);
    u_cosA_val = p(9);
    u_cosO_val = p(10);
 
    ut    = exp(-t/u_expT_val)*(u_sinA_val*sin(u_sinO_val*t)+u_cosA_val*cos(u_cosO_val*t));
    dxdt = [ x(2); ...
            -alpha1_val*x(2)-alpha0_val*x(1)*(1+gamma1_val*(x(1)^2))+...
                beta0_val*(1+gamma2_val*x(1))*ut ];
 
end % function

The Third Order Nonlinear Difference Equation: the Difference Equation as a MATLAB Function

This model can be expressed in state-space form as

where

with initial conditions

and the input sequence

function [y_out, x_out] = Solve_ORD3_NONLIN_DifferenceEquation(timeT, x0, u_in )
    persistent x1k x2k x3k uk_1
 
    % --- let's initialize the persistent variables ---
    if isempty(x1k)|isempty(x2k) | isempty(x3k)
        x1k = x0(1);
        x2k = x0(2);
        x3k = x0(3);
    end
    if isempty(uk_1)
        uk_1 = 0;
    end
     % --- let's initialize the persistent variables ---
    
      nT = numel(timeT); % how many time steps?
 
     % --- let's initialize the output arrays ---
     y_out = zeros(nT,1);
     x_out = zeros(nT, 3);
     % --- let's initialize the output arrays ---
 
     % --- the simulation loop ---
     for k =1:nT
        uk    = u_in(k);
        x1kp1 = x2k;
        x2kp1 = x3k;
        x3kp1 = (x3k*x2k*x1k*(x1k-1)*uk_1+uk)/(1+x2k^2+x1k^2);
        yk    = x3k;
 
        % store the results
        y_out(k) = yk;
        x_out(k,:) = [x1k; x2k; x3k];
 
        % update state and past-input variables (the persistent variables)
        x1k = x1kp1;
        x2k = x2kp1;
        x3k = x3kp1;
        uk_1 = uk;
 
     end % for k
     % --- the simulation loop ---
 
end % function