ctrb function returns NaN or inf when dealing with large systems
Show older comments
I’m currently working with continuous-time, large state-space systems — specifically, systems with a number of states on the order of 10^3, and numbers of inputs and outputs on the order of 10^2.
I’m trying to analyze the reachability and observability properties of these systems. To that end, I’ve been using MATLAB’s ctrb function, which, according to the documentation, returns the reachability matrix. However, I’ve encountered a major issue: the resulting matrix (Co) contains a significant number of NaN and Inf entries.
For reference, in one example the matrix Co has dimensions 1,624 x 311,808, resulting in a total of 506,376,192 elements. Of these:
- 466,771,919 entries (92.18%) are NaN, counted with sum(isnan(Co), 'all')
- 37,732 entries (0.0075%) are Inf or -Inf, counted with sum(isinf(Co), 'all')
Has anyone encountered a similar issue when working with large-scale systems?
How did you handle it or work around it?
Any suggestions, tips, or alternative methods would be highly appreciated!
9 Comments
Just a test. It seems that the number of states cannot exceed 348!
Case 1:
n = 348; % number of states
p = - ones(1, n); % 100% stable poles
G = zpk([], p, prod(p)); % transfer function
% sys = canon(G, 'companion') % The "canon" command failed because the controllability matrix is nearly singular.
sys = compreal(G); % convert ot state-space
rk = rank(ctrb(sys.A, sys.B))
Case 2:
n = 349; % number of states
p = - ones(1, n); % 100% stable poles
G = zpk([], p, prod(p)); % transfer function
sys = compreal(G); % convert ot state-space
rk = rank(ctrb(sys.A, sys.B))
Case 3: An unexplainable result?
n = 250; % number of states
sys = rss(n); % generate random continuous-time state-space
Tf = isstable(sys)
rk = rank(ctrb(sys.A, sys.B))
Mathieu NOE
on 7 Apr 2025
Basing on the idea of @Sam Chak, I run the following script multiple time (also with different n). It seems like (but I am not sure at all) that it is correlated to the fact that the matrix A (and B if more than 1 input) is bad conditioned and 
for a certain 
degenerates to values close to infinity. I don't have any idea on how to treat that, a part from checking manually the presence of numerical schemes that leads to, at least, two rows linearly dependent and hence concluding the matrix Co can never be full rank. However this is unfeasible with problems like mine, where size(B) is ~10^3 x 10^2.
for a certain degenerates to values close to infinity. I don't have any idea on how to treat that, a part from checking manually the presence of numerical schemes that leads to, at least, two rows linearly dependent and hence concluding the matrix Co can never be full rank. However this is unfeasible with problems like mine, where size(B) is ~10^3 x 10^2.n = 250;
sys = rss(n);
Tf = isstable(sys)
try
rk = rank(ctrb(sys.A, sys.B)) % Check reachability with Matlab function
catch ME
warning("ctrb gave the following error: \n\t %s \n", ME.message)
end
% Check reachability by hand
A = sys.A; B = sys.B; % Scomposing them for simplicity
Aprod = sys.A;
Co = [B, Aprod*B];
for ii = 2:n-1
fprintf("n = %d, cond(Aprod) = %.3e\n", ii-1, cond(Aprod));
Aprod = Aprod*A; % A^n = A^(n-1)*A
product = Aprod*B; % in C0 goes A^n*B
try
assert(sum(isnan(product), 'all')==0, 'Elements in product are NaN: #: %d, iteration: %d, rank(C0) = %d', ...
sum(isnan(product), 'all'), ii, rank(Co)); % Check if NaN
assert(sum(isinf(product), 'all')==0, 'Elements in product are Inf: #: %d, iteration: %d, rank(C0) = %d', ...
sum(isinf(product), 'all'), ii, rank(Co)); % Check if Inf
catch e
disp(product) % Show where the cycle broke
rethrow(e)
end
Co = [Co, product];
if rank(Co) == n, fprintf('Iterations: %d', ii); break; end
end
rk = rank(Co)
Corby
on 7 Apr 2025
Mathieu NOE
on 7 Apr 2025
Sure , a reduced order system cannot be compared or taken as 100% "equivalent" to the full system. Now the engineer must ask himself what does it mean to remove the poorly reachable / observable states ? Is it still physically meaningfull ?
Case 5: The ctrb() function can handle systems with a number of states greater than 1,000.
n = 2000; % number of states
A0 = zeros(n, 1);
A1 = eye(n);
A1(:,n) = [];
A = [A0, A1]; % State matrix
B = zeros(sqrt(numel(A)), 1);
B(end) = 1; % Input matrix
Co = ctrb(A, B); % Controllability matrix
rk = rank(Co) % expect to return rk = n
Hi Corby,
I think the fundamental problem that forming the controllability matrix necessarily means raising the system matrix to a large power, which at some point will overflow if the spectral radius is greater than unity and n is large enough. I think you touched on this issue in this comment, though I don't think that necessarily means the system matrix is poorly conditioned.
For example
rng(100);
n = 250;
sys = rss(n,1,5);
size(sys)
rss should generate a stable system
Tf = isstable(sys)
but it doesn't because
max(real(eig(sys.a)))
So shift the A matrix to ensure stability (do you know your system is stable?)
sys.a = sys.a - eye(n);
max(real(eig(sys.a)))
We see that raising the system matrix to the i'th power breaks down at i = 195
for ii = 1:n
if any(isnan(sys.A^ii),'all') || any(isinf(sys.A^ii),'all')
break
end
end
ii
At i = 194, elements are getting close to realmax
[max(sys.A^(ii-1),[],'all') , min(sys.A^(ii-1),[],'all'), realmax]
The spectral radius of the system matrix is
r = max(abs(eig(sys.A)))
The spectral radius of A^n is r^n, so we might estimate the maximum achievable n as
nlimit = log(realmax)/log(r)
which at least tracks nicely with the result above, though I don't really know how rigorous that really is.
Anyway, that leaves us searching for other ways to test for controllability.
One option is to compute the controllability Gramian
Wc = gram(sys,'c');
which at least is returning a finite result
all(isfinite(Wc),'all')
By definiton the Gramian is symmetric
issymmetric(Wc)
It it positive definite if (A,B) are a controllable pair.
The simple test would be
e = eig(Wc);
all(imag(e)==0)
min(e)
Taken literally, that means the systme is not controllable (would that be miraculous for a randomly generated system?).
try
C = chol(Wc)
catch ME
ME.message
end
Does that mean they system isn't controllable? Or does that mean we are suffering from numerical issues in one of the steps along the way?
Another option might be to tranform the system to modal form. I think a sufficient condition for controllability is that every row of the transformed B matrix has at least one non-zero entry (check that)
sysm = modalreal(sys);
all(any(sysm.b,2))
Of course, here we could have the problem where a value is very close to zero. Here that's not the case.
min(min(abs(sysm.b),[],2))
In any case, there might be options besides checking the controllability matrix. Also, consider looking into other numerical conditioning schemes.
Seems like a hard problem in general.
If you don't mind, what kind of system is being modeled?
More importantly, why is it important to know if the system is controllable? What would that information be used for if it could be accurately ascertained?
Corby
on 8 Apr 2025
Accepted Answer
Hi @Corby
Thank you and Paul for explaining why the ctrb() function fails for very high-order LTI systems when the spectral radius of the state matrix 
exceeds the realmax value. The Controllability Gramian is also an effective approach to assess the system's controllability through the Lyapunov equation.
exceeds the realmax value. The Controllability Gramian is also an effective approach to assess the system's controllability through the Lyapunov equation. However, when dealing with very high-order LTI systems, we can use the Hautus Lemma, commonly known as the Popov-Belevitch-Hautus test, not only to determine the controllability of the system but also to identify the minimum number of control inputs required to render the LTI system controllable. I highly recommend watching Prof. Steve Brunton's video.
Simple code for the Popov-Belevitch-Hautus (PBH) test:
%% State-space system
n = 500; % number of system states
sys = rss(n); % random state-space
%% Setting for PBH Test
tol = 1e-16; % tolerance to use in the rank computation
%% Call isctrb() to determine controllability
TF = isctrb(sys, tol)
%% Code for Popov-Belevitch-Hautus (PBH) test
function TF = isctrb(sys, tol)
% Extract state and input matrices from the system
A = sys.A; % state matrix from sys
B = sys.B; % input matrix from sys
% Find eigenvalues of the system
p = eig(A); % eigenvalues (poles)
ns = sqrt(numel(A)); % number of states
I = eye(ns); % identity matrix
% Initialize rank vector
rk = zeros(ns, 1);
% Perform PBH test
for i = 1:ns
rk(i) = rank([A-p(i)*I, B], tol);
% can break here; not need to loop all
if rk(i) ~= ns
break
end
end
% Display result
if any(rk ~= ns, 'all')
TF = 0;
disp('The system is not controllable.');
else
TF = 1;
disp('The system is controllable.');
end
end
2 Comments
Paul
on 8 Apr 2025
Good idea to investigate the PBH test.
Would there be any concerns with determining rank of the test matrix, which involves a tolerance? Because we don't need to know the exact rank, only whether or not it's full rank, perhaps it would be better to just look at its minimum singular value (svd), but that will still involve a judgement against a tolerance. I have no idea what to expect for accuracy of svd for matrices of this scale, especially for the minimum singular value.
Sam Chak
on 9 Apr 2025
Hi @Paul, thank you for your advice. The code was written hastily yesterday. I have tested it, and the rank computation indeed has issues for systems with ill-conditioned state matrices (after being transformed from equivalent well-conditioned matrices). Allowing users to specify the tolerance helps to alleviate this issue.
Both of you are welcome to improve the code as needed.
More Answers (0)
Categories
Find more on Matrix Computations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
