ctrb function returns NaN or inf when dealing with large systems

25 visualizaciones (últimos 30 días)
Corby
Corby el 7 de Abr. de 2025
Comentada: Sam Chak el 9 de Abr. de 2025
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 comentarios
Mostrar 7 comentarios más antiguos
Paul
Paul el 7 de Abr. de 2025
Editada: Paul el 8 de Abr. de 2025
Ran in:
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)
State-space model with 1 outputs, 5 inputs, and 250 states.
rss should generate a stable system
Tf = isstable(sys)
Tf = logical
0
but it doesn't because
max(real(eig(sys.a)))
ans = 7.6458e-14
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)))
ans = -1.0000
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
ii = 195
At i = 194, elements are getting close to realmax
[max(sys.A^(ii-1),[],'all') , min(sys.A^(ii-1),[],'all'), realmax]
ans = 1×3
1.0e+308 * 0.1267 -0.1341 1.7977
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The spectral radius of the system matrix is
r = max(abs(eig(sys.A)))
r = 38.4665
The spectral radius of A^n is r^n, so we might estimate the maximum achievable n as
nlimit = log(realmax)/log(r)
nlimit = 194.4723
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')
ans = logical
1
By definiton the Gramian is symmetric
issymmetric(Wc)
ans = logical
1
It it positive definite if (A,B) are a controllable pair.
The simple test would be
e = eig(Wc);
all(imag(e)==0)
ans = logical
1
min(e)
ans = -6.4753e-13
Taken literally, that means the systme is not controllable (would that be miraculous for a randomly generated system?).
Or we can use chol to test, and it too fails.
try
C = chol(Wc)
catch ME
ME.message
end
ans = 'Matrix must be positive definite.'
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))
ans = logical
1
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))
ans = 0.0122
Or try a balanced realization (balreal)
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
Corby el 8 de Abr. de 2025
Dear @Paul, Thank you for your kind answer. The Gramian seems like a very clever way to assess controllability.
The system in question models the eddy current response of a 3D structure to external magnetic field variations.
Controllability is of interest because, if confirmed, it would imply that a desired eddy current pattern can be achieved by controlling a finite set of currents.

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuesta aceptada

Sam Chak
Sam Chak el 8 de Abr. de 2025
Editada: Sam Chak el 9 de Abr. de 2025
Ran in:
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.
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)
The system is controllable.
TF = 1
%% 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 comentarios
Paul
Paul el 8 de Abr. de 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
Sam Chak el 9 de Abr. de 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.
Hi @Corby, I have updated the code for the PBH test. Please take a look if you are interested.
Both of you are welcome to improve the code as needed.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Matrix Computations en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2024b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by