%[text] %[text:anchor:T_495B1013] # Equilibrium States of LTI Dynamic Systems
%[text] 
%[text:tableOfContents]{"heading":"Table of Contents"}
%[text] %[text:anchor:H_0F439577] ## Introduction
%[text] 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.
%[text] 
%[text] What you will learn: 
%[text] -  How to implement and solve the equilibrium state equation $(I-A) \\bar x = B \\bar u$ using Matlab.
%[text] - How to cope with the critical cases, when no equilibrium or infinite equilibrium states exist. \
%[text] 
%[text] %[text:anchor:H_333592A2] ## Equilibrium State Equation of LTI Systems
%[text] Consider a LTI dynamic system:
%[text]{"align":"center"} $\\left \\{\n\\begin{array}{l}\nx(k+1) =\nA x(k) + B u(k)\\\\\n\\\\\ny(k) = C x(k) + D u(k)\n\\end{array}\n\\right . $
%[text] and consider a constant input sequence 
%[text]{"align":"center"} $u(k) = \\bar u , \\quad k \\ge 0$
%[text] 
%[text] Therefore, the following equation must be solved to find the equilibrium state $\\bar x$, if any, associated with the constant input $\\bar u$`:`
%[text]{"align":"center"} $\\bar x = A \\bar x + B \\bar u \\, \\Longrightarrow \\,\n(I-A) \\bar x = B \\bar u$
%[text] The following two cases have to be considered:
%[text] - ${\\rm det} \\, (I-A) \\not= 0$
%[text] - ${\\rm det} \\, (I-A) = 0$ \
%[text] How is this equation solved using Matlab?
%[text] 
%[text] %[text:anchor:H_D9726AE7] ### Notation
%[text] Given an $n\\,\\text{x}\\,n$ matrix $&dollar&;M&dollar&;$,  define
%[text] - the **column space** of $M\n$is the span of its column vectors: $M = \\left\[ m\_1 \\, m\_2\\,\\cdots\\, m\_n \\right\]\\,,\\; m\_i \\in \\mathbb{R}^{n} \\;i=1,\\,2,\\ldots\\, n\\; \\Rightarrow \\text{colsp} \\left(M\\right) = \\langle m\_1 \\, m\_2\\,\\cdots\\, m\_n \\rangle$
%[text] - the **null space** or kernel of $M$is the subspace of $\\mathbb{R}^{n}$ which is mapped to the zero vector: $\\text{Null}\\left(M\\right) = \\text{ker}\\left(M\\right) = \\left\\{x\\in \\mathbb{R}^{n}\\,:\\; Ax=0\\right\\}$ \
%%
%[text] %[text:anchor:H_9D09ECFD] ## Solving the Equilibrium State Equation: One and Only One Solution Exists 
%[text] %[text:anchor:H_979900BF] We are dealing with the case when ${\\rm det} \\, (I-A) \\not= 0$
%[text] Cleaning up the Matlab workspace
close all
clear
clc

%[text] 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) %[output:21b7e72b]
%%
%[text]  Show the eigenvalues of the matrix$A$ 
lambdaA = eig(LTID_sys.A) %[output:0734ad8f]
%[text] Check if ${\\rm det} \\, (I-A) \\not= 0$
assert(rank(eye(SysOrder)-LTID_sys.A)==SysOrder)
% if the logical condition is not met, an error is thrown
%%
%[text] Now, the unique equilibrium state $\\bar x$ is 
%[text]{"align":"center"} $\\bar x = \\left(I-A\\right)^{-1}\\,B\\,\\bar u$
% 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 %[output:07ea9d70]
% 2nd approach - using the linsolve function 
%[text] Let's solve the linear algebra problem
%[text]  $(I-A) \\bar x = B \\bar u \\; \\Longleftrightarrow \\; M\\,v = b$
%[text] exploiting the MATLAB function **`linsolve`**
%[text] ```matlabCodeExample
%[text] vSOL = linsolve(M, b)
%[text] ```
%[text] 
% Avoid the matrix inversion: numerically more stable and efficient
% approach
barX_bis = linsolve(eye(SysOrder)-LTID_sys.A, LTID_sys.B*barU) %[output:76139499]
%[text] Obviously, it's the same solution. The second approach is more reliable, so it's preferable.
%[text] Given an equilibrium state $\\bar x$ and the corresponding equilibrium input $\\bar u$, the equilibrium output is given by: 
%[text]{"align":"center"} $\\bar y = \\left\[ C \\left(I-A\\right)^{-1}\\,B+D \\right\] \\, \\bar u$
%[text] **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 %[output:1a47cb82]
%[text] **Note**: an alternative, more numerically reliable, expression for the computation of the equilibrium output $\\bar y$ is simply
%[text]{"align":"center"} $\\bar y = C \\, \\bar x + D \\, \\bar u$
barY_ALT = LTID_sys.C * barX_bis + LTID_sys.D * barU %[output:8258e2c7]
%%
%[text] %[text:anchor:H_AE69A127] ## Solving the Equilibrium State Equation: No Solution or Infinite Solutions 
%[text] 
%[text] %[text:anchor:H_3B13BE78] We are dealing with the case ${\\rm det} \\, (I-A) = 0$
%[text] Cleaning up the Matlab workspace
close all
clear
clc

SysOrder = 6; % the system order 
%%
%[text] Consider a generic SISO discrete-time LTI system with **at least** a simple eigenvalue in $z=+1$
foundIT = false;
i=0;
while ~foundIT %[output:group:0b8d261c]
    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! %[output:8a691dc1]
    else
        fprintf(1,'trial n. %d\n', i); %[output:9d9072ef]
    end

end %[output:group:0b8d261c]
%%
% 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
%[text] The following two situations may occur:
%[text] 1. $B\\,\\bar u \\in \\text{colsp}\\left(I-A \\right)$ : $&dollar&;\\exists \\, \\infty \\, \\, \\textrm{equilibrium states} \\, \\, \\bar x&dollar&;, &dollar&;\\exists \\, \\infty \\, \\, \\textrm{equilibrium outputs} \\, \\, \\bar y&dollar&;$
%[text] 2. $B\\,\\bar u \\notin \\text{colsp}\\left( I-A \\right)$ : $&dollar&;\\not\\exists \\, \\, \\textrm{equilibrium states} \\, \\, \\bar x&dollar&;, &dollar&;\\not\\exists \\, \\, \\textrm{equilibrium outputs} \\, \\, \\bar y&dollar&;$ \
%%
%[text] 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.
%[text] - 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 <<---

%[text] - 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. \
%%
%[text] %[text:anchor:H_665277FB] ###  Check the Existence of Solutions
%[text] %[text:anchor:H_C4710095] Obviously, a solution exists if $\\text{colsp} \\left(B \\right) \\in \\text{colsp}\\left( I-A \\right)$
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) %[output:group:7d069409]
    % 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!') %[output:8f4ccca0]
    
    % solving the problem using linsolve
    X0 = linsolve(ImA, bb) %[output:5fa88416] %[output:385f7fc6]
    % evaluate the accuracy of the found solution
    SolAccuracy = norm(ImA*X0-bb) %[output:50099f80]
    % 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 %[output:4577831c]

    sympref(old_prefs);
end % if (check_rho_IAB > rhoIA) %[output:group:7d069409]
%[text] 
%[text] 
%[text] 

%[appendix]{"version":"1.0"}
%---
%[metadata:view]
%   data: {"layout":"inline","rightPanelPercent":40}
%---
%[output:21b7e72b]
%   data: {"dataType":"textualVariable","outputData":{"name":"barU","value":"0.8351"}}
%---
%[output:0734ad8f]
%   data: {"dataType":"matrix","outputData":{"columns":1,"name":"lambdaA","rows":6,"type":"complex","value":[["0.8680 + 0.0000i"],["0.6983 + 0.0000i"],["0.3115 + 0.0000i"],["-0.7858 + 0.1005i"],["-0.7858 - 0.1005i"],["-0.9286 + 0.0000i"]]}}
%---
%[output:07ea9d70]
%   data: {"dataType":"matrix","outputData":{"columns":1,"name":"bar_X","rows":6,"type":"double","value":[["-9.4920"],["-2.0792"],["4.0460"],["-7.3828"],["7.9059"],["3.1190"]]}}
%---
%[output:76139499]
%   data: {"dataType":"matrix","outputData":{"columns":1,"name":"barX_bis","rows":6,"type":"double","value":[["-9.4920"],["-2.0792"],["4.0460"],["-7.3828"],["7.9059"],["3.1190"]]}}
%---
%[output:1a47cb82]
%   data: {"dataType":"textualVariable","outputData":{"name":"barY","value":"-5.1567"}}
%---
%[output:8258e2c7]
%   data: {"dataType":"textualVariable","outputData":{"name":"barY_ALT","value":"-5.1567"}}
%---
%[output:9d9072ef]
%   data: {"dataType":"text","outputData":{"text":"trial n. 1\ntrial n. 2\ntrial n. 3\ntrial n. 4\ntrial n. 5\ntrial n. 6\ntrial n. 7\ntrial n. 8\ntrial n. 9\ntrial n. 10\ntrial n. 11\ntrial n. 12\ntrial n. 13\ntrial n. 14\ntrial n. 15\ntrial n. 16\ntrial n. 17\ntrial n. 18\ntrial n. 19\ntrial n. 20\ntrial n. 21\ntrial n. 22\ntrial n. 23\ntrial n. 24\ntrial n. 25\ntrial n. 26\ntrial n. 27\ntrial n. 28\ntrial n. 29\ntrial n. 30\ntrial n. 31\ntrial n. 32\ntrial n. 33\ntrial n. 34\ntrial n. 35\ntrial n. 36\ntrial n. 37\ntrial n. 38\ntrial n. 39\ntrial n. 40\ntrial n. 41\ntrial n. 42\ntrial n. 43\ntrial n. 44\ntrial n. 45\ntrial n. 46\ntrial n. 47\ntrial n. 48\ntrial n. 49\ntrial n. 50\ntrial n. 51\ntrial n. 52\ntrial n. 53\ntrial n. 54\ntrial n. 55\ntrial n. 56\ntrial n. 57\ntrial n. 58\ntrial n. 59\ntrial n. 60\ntrial n. 61\ntrial n. 62\ntrial n. 63\ntrial n. 64\ntrial n. 65\ntrial n. 66\ntrial n. 67\ntrial n. 68\ntrial n. 69\ntrial n. 70\ntrial n. 71\ntrial n. 72\ntrial n. 73\ntrial n. 74\ntrial n. 75\ntrial n. 76\ntrial n. 77\ntrial n. 78\ntrial n. 79\ntrial n. 80\ntrial n. 81\ntrial n. 82\ntrial n. 83\ntrial n. 84\ntrial n. 85\ntrial n. 86\ntrial n. 87\ntrial n. 88\ntrial n. 89\ntrial n. 90\ntrial n. 91\ntrial n. 92\ntrial n. 93\ntrial n. 94\ntrial n. 95\ntrial n. 96\ntrial n. 97\ntrial n. 98\ntrial n. 99\ntrial n. 100\ntrial n. 101\ntrial n. 102\ntrial n. 103\ntrial n. 104\ntrial n. 105\ntrial n. 106\ntrial n. 107\ntrial n. 108\ntrial n. 109\ntrial n. 110\ntrial n. 111\ntrial n. 112\ntrial n. 113\ntrial n. 114\ntrial n. 115\ntrial n. 116\ntrial n. 117\ntrial n. 118\ntrial n. 119\ntrial n. 120\ntrial n. 121\ntrial n. 122\ntrial n. 123\ntrial n. 124\ntrial n. 125\ntrial n. 126\ntrial n. 127\ntrial n. 128\ntrial n. 129\ntrial n. 130\ntrial n. 131\ntrial n. 132\ntrial n. 133\ntrial n. 134\ntrial n. 135\ntrial n. 136\ntrial n. 137\ntrial n. 138\ntrial n. 139\ntrial n. 140\ntrial n. 141\ntrial n. 142\ntrial n. 143\ntrial n. 144\ntrial n. 145\ntrial n. 146\ntrial n. 147\ntrial n. 148\ntrial n. 149\ntrial n. 150\ntrial n. 151\ntrial n. 152\ntrial n. 153\ntrial n. 154\ntrial n. 155\ntrial n. 156\ntrial n. 157\ntrial n. 158\ntrial n. 159\ntrial n. 160\ntrial n. 161\ntrial n. 162\ntrial n. 163\ntrial n. 164\ntrial n. 165\ntrial n. 166\ntrial n. 167\ntrial n. 168\ntrial n. 169\ntrial n. 170\ntrial n. 171\ntrial n. 172\ntrial n. 173\ntrial n. 174\ntrial n. 175\ntrial n. 176\ntrial n. 177\ntrial n. 178\ntrial n. 179\ntrial n. 180\ntrial n. 181\ntrial n. 182\ntrial n. 183\ntrial n. 184\ntrial n. 185\ntrial n. 186\ntrial n. 187\ntrial n. 188\ntrial n. 189\ntrial n. 190\ntrial n. 191\ntrial n. 192\ntrial n. 193\ntrial n. 194\ntrial n. 195\ntrial n. 196\ntrial n. 197\ntrial n. 198\ntrial n. 199\ntrial n. 200\ntrial n. 201\ntrial n. 202\ntrial n. 203\ntrial n. 204\ntrial n. 205\ntrial n. 206\ntrial n. 207\ntrial n. 208\ntrial n. 209\ntrial n. 210\ntrial n. 211\ntrial n. 212\ntrial n. 213\ntrial n. 214\ntrial n. 215\ntrial n. 216\ntrial n. 217\ntrial n. 218\ntrial n. 219\ntrial n. 220\ntrial n. 221\n","truncated":false}}
%---
%[output:8a691dc1]
%   data: {"dataType":"text","outputData":{"text":"Eureka!\n","truncated":false}}
%---
%[output:8f4ccca0]
%   data: {"dataType":"text","outputData":{"text":"Infinite equilibrium states and outputs!\n","truncated":false}}
%---
%[output:5fa88416]
%   data: {"dataType":"warning","outputData":{"text":"Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.131522e-17."}}
%---
%[output:385f7fc6]
%   data: {"dataType":"matrix","outputData":{"columns":1,"name":"X0","rows":6,"type":"double","value":[["-1.0718"],["-0.3311"],["-0.2159"],["-0.0228"],["0.4825"],["0"]]}}
%---
%[output:50099f80]
%   data: {"dataType":"textualVariable","outputData":{"name":"SolAccuracy","value":"1.6329e-16"}}
%---
%[output:4577831c]
%   data: {"dataType":"symbolic","outputData":{"name":"X_SET","value":"\\left\\lbrack \\begin{array}{c}\n0.7367\\,p_1 -1.0718\\\\\n-0.0534\\,p_1 -0.3311\\\\\n0.1966\\,p_1 -0.2159\\\\\n-0.0547\\,p_1 -0.0228\\\\\n0.4825-0.3267\\,p_1 \\\\\n0.5532\\,p_1 \n\\end{array}\\right\\rbrack"}}
%---