barU = 0.8351
Equilibrium States of LTI Dynamic Systems
Table of Contents
Introduction
In this live script, we illustrate how to use Matlab to determine the equilibrium state of an LTI system, if one exists, given a constant input. We also show how to handle cases where no equilibrium or infinite equilibrium states exist.
What you will learn:
- How to implement and solve the equilibrium state equation
using Matlab. - How to cope with the critical cases, when no equilibrium or infinite equilibrium states exist.
Equilibrium State Equation of LTI Systems
Consider a LTI dynamic system:
and consider a constant input sequence
Therefore, the following equation must be solved to find the equilibrium state 
, if any, associated with the constant input 
:
, if any, associated with the constant input :
The following two cases have to be considered:
How is this equation solved using Matlab?
Notation
Given an 
matrix M, define
matrix M, define- the column space of Mis the span of its column vectors:
- the null space or kernel of Mis the subspace of
which is mapped to the zero vector: 
Solving the Equilibrium State Equation: One and Only One Solution Exists
We are dealing with the case when 
Cleaning up the Matlab workspace
close all
clear
clc
Create a generic SISO discrete-time LTI system, using the state-space representation
SysOrder = 6; % the system order
rng(0); % configuring the random number generator
LTID_sys = drss(SysOrder); % a SISO discrete time system of order SysOrder
barU = randn(1) % select a random equilibrium input (constant input)
Show the eigenvalues of the matrixA
lambdaA = eig(LTID_sys.A)
Check if
assert(rank(eye(SysOrder)-LTID_sys.A)==SysOrder)
% if the logical condition is not met, an error is thrown
Now, the unique equilibrium state is
is 
% solve the equation - 1st approach: solution by using the inverse matrix
% <-- this could be an inefficient approach, from the numerical point of view !!
% Indeed, the matrix inversion algorithm could be ill conditioned: refer,
% for example, to
%
% https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010/resources/lecture-3-multiplication-and-inverse-matrices/
%
bar_X = inv(eye(SysOrder)-LTID_sys.A)*LTID_sys.B*barU
% 2nd approach - using the linsolve function
Let's solve the linear algebra problem
exploiting the MATLAB function linsolve
vSOL = linsolve(M, b)
% Avoid the matrix inversion: numerically more stable and efficient
% approach
barX_bis = linsolve(eye(SysOrder)-LTID_sys.A, LTID_sys.B*barU)
Obviously, it's the same solution. The second approach is more reliable, so it's preferable.
Given an equilibrium state and the corresponding equilibrium input , the equilibrium output is given by:
and the corresponding equilibrium input , the equilibrium output is given by: 
Remark: the formula relies on the inversion of a matrix. Thus, once more, the result could be ill-conditioned.
barY = (LTID_sys.C*inv(eye(SysOrder)-LTID_sys.A)*LTID_sys.B+LTID_sys.D)*barU
Note: an alternative, more numerically reliable, expression for the computation of the equilibrium output 
is simply
is simply
barY_ALT = LTID_sys.C * barX_bis + LTID_sys.D * barU
Solving the Equilibrium State Equation: No Solution or Infinite Solutions
We are dealing with the case 
Cleaning up the Matlab workspace
close all
clear
clc
SysOrder = 6; % the system order
Consider a generic SISO discrete-time LTI system with at least a simple eigenvalue in 
foundIT = false;
i=0;
while ~foundIT
i=i+1;
dummySYS = drss(SysOrder);
lambdaA = eig(dummySYS.A);
if any(lambdaA==+1) % <-- look for an eigenvalue in z=+1
foundIT = true;
disp('Eureka!'); % found it!
else
fprintf(1,'trial n. %d\n', i);
end
end
% store the LTI system in LTID_sys_eig1
LTID_sys_eig1 = dummySYS;
% now select a constant input
barU = randn(size(LTID_sys_eig1.B,2)); % # of columns in B is equal to # of inputs
The following two situations may occur:
: 
: 
The actual code (with commented the piece of code starting from the 50th to the 56th line) copes with the second case: what if there is no equilibrium at all.
- Leave commented out the following piece of code and execute the code in the following section: what happens?
% uncomment the following piece of code,
% in order to have infinite equilibrium states
% the code forces colsp(B) to be a subspace of colsp(I-A)
% -->> uncomment starting from the line below <<---
A = LTID_sys_eig1.A;
K = randn(SysOrder,1); % select n real coeffs, where n is #
% of column vectors in the matrix A
B = (eye(size(A))-A)*K; % Now B is a linear combination of the column vectors in A
% NB B is a column vector - the LTID_sys_eig1 system is a SISO
% LTI system by construction
LTID_sys_eig1.B = B; % forcing the B matrix of the LTI system
% -->> uncomment the piece of code above this line <<---
- Then, uncomment the code, execute it and run again the code in the following section: this time, the code copes with the first situation: there are infinite solutions to the equilibrium problem. Have a look on how to find the whole set of solutions using MATLAB.
Check the Existence of Solutions
Obviously, a solution exists if 
matA = LTID_sys_eig1.A;
matB = LTID_sys_eig1.B;
bb = matB*barU;
ImA = eye(size(matA))-matA; % (I-A)
rhoIA = rank(ImA);
rhoB = rank(matB);
check_rho_IAB = rank([ImA, matB]);
if (check_rho_IAB > rhoIA)
% the column space of B is NOT in the column space by (I-A)
disp('There is no equilibrium state and no equilibrium output!')
% If you try to solve the linear matrix equation anyway,
% you receive a nonsense as result and a warning message
% -- just uncomment the following line to have a look
nonsenseXsolution = linsolve(ImA, bb)
else
% the column space of B IS in the column space by (I-A)
disp('Infinite equilibrium states and outputs!')
% solving the problem using linsolve
X0 = linsolve(ImA, bb)
% evaluate the accuracy of the found solution
SolAccuracy = norm(ImA*X0-bb)
% Note: a very small error norm
% Let's write the set of solutions as linear combination of the
% particular found solution and the basis of the null space Null(I-A)
null_I_A = null(ImA); % a basis of the null space Null(I-A)
nP = size(null_I_A,2); % how many column vectors?
% the whole set of equilibrium states as:
% general solution = Xeq_linS + freeP * nV
%
% nV in Null(I-A), p is a set of free parameters, as many as the
% dimension of the null space Null(I-A)
old_prefs = sympref(); % store the actual Symbolic Toolbox preferences
sympref('MatrixWithSquareBrackets',true); % matrix style
sympref('FloatingPointOutput',true); % To display the results in floating-point format
syms p [1 nP]
X_SET = X0 + null_I_A*p
sympref(old_prefs);
end % if (check_rho_IAB > rhoIA)
