[Math] Minimal realization of a state space model

control theorylinear-controlMATLABoptimal control

Minimal realization of a transfer function is about cancelling out poles against zeros.

But the minimal realization of a state space model is about cancelling out non-controllable and non-observable states.

My question is how I can do that? Can I use PBH (Popov, Belevich, Hautus)-test to find which each eigenvalue gives non-controllable or observable state?

$$\operatorname{rank} ([(A-\lambda I) B]) = n \forall \lambda $$

$$\operatorname{rank} ([(A-\lambda I); C]) = n \forall \lambda $$

After I have found the states who are uncontrollable or unobservable. Can I then create a new state space model without those states?

Best Answer

The easiest way of finding the minimal state space model I think would be using the Kalman decomposition. This allows you to find a similarity transformation that makes it easy to split the state space model into a minimal (controllable and observable) and non-minimal (not controllable or not observable) form.

When calculating this decomposition you could use the Hautus test, but you could also use the controllability and observability matrix. I made a quick implementation of both in MATLAB, whose code can be seen below.

Hautus test:

function [M, Ak, Bk, Ck, N] = kalman_decomp_hautus(A, B, C)

[V,D] = eig(A);
n = length(A);

M1 = []; M2 = []; M3 = []; M4 = [];
N = zeros(1, 4);
for k = 1 : n
    tempC = rank([A - D(k,k) * eye(n) B]) - n;
    tempO = rank([A - D(k,k) * eye(n); C]) - n;
    if tempC == 0
        if tempO == 0
            M2 = [M2 V(:,k)];
            N(2) = N(2) + 1;
        else
            M1 = [M1 V(:,k)];
            N(1) = N(1) + 1;
        end
    else
        if tempO == 0
            M4 = [M4 V(:,k)];
            N(4) = N(4) + 1;
        else
            M3 = [M3 V(:,k)];
            N(3) = N(3) + 1;
        end
    end
end
M = [M1 M2 M3 M4];

Ak = M \ A * M;
Bk = M \ B;
Ck = C * M;

Using this implementation does have the downside that the decomposed state space model matrices can have complex numbers, because if the original $A$ matrix (assuming it is real valued) has complex conjugate eigenvalues then its eigenvectors will be complex conjugate as well, and therefore the similarity transformation also. You could of course try to detect this, but I just wanted a working example.

Controllability and observability matrix:

function [M, Ak, Bk, Ck, N] = kalman_decomp_matrix(A, B, C, tol)

cc = ctrb(A, B);
oo = obsv(A, C);

Sc = rref(cc')'; Sc = Sc(:,1:rank(cc));
Nc = null(cc');
So = rref(oo)'; So = So(:,1:rank(oo));
No = null(oo);

M1 = rref(round(Sc * (Sc \ No),tol)')';
N(1) = rank(M1);
M1 = M1(:,1:N(1));
M2 = rref(round(Sc * (Sc \ So),tol)')';
N(2) = rank(M2);
M2 = M2(:,1:N(2));
M3 = rref(round(Nc * (Nc \ No),tol)')';
N(3) = rank(M3);
M3 = M3(:,1:N(3));
M4 = rref(round(Nc * (Nc \ So),tol)')';
N(4) = rank(M4);
M4 = M4(:,1:N(4)); 
M = [M1 M2 M3 M4];

Ak = M \ A * M;
Bk = M \ B;
Ck = C * M;

I have not done many tests, but with one test I noticed that finding the intersect of two spans of vectors could give wrong results due to limited numerical accuracy, so I added rounding (using tol=13 seems to omit errors in my case).

Once you have your Kalman decomposition then the minimal realisation of the state space model can be constructed using:

list = N(1) + 1:N(1) + N(2);
Am = Ak(list,list);
Bm = Bk(list,:);
Cm = Ck(:,list);
Dm = D;

If you just want to minimal realisation then you might be able to calculate it faster. But you do need to ensure that the total $M$ matrix is square and full rank.

Related Question