I'm trying to program the z-transform in wxMaxima which doesn't have it programmed but not by definition but by using the Scilab approach. Scilab to calculate the z-transform first converts the transfer function to the state space, after that the system must be discretized and after that converted to z transfer function, I need this because of some algebraic calculations that I need to do to analyze stability of a system in function of the sample period.
Right now I'm stranded with the function balanc() which finds a similarity transform such that
Ab = X^(-1) . A . X
as approximately equal row and column norms.
Most of my code in wxMaxima to reach in the near future has been done by translating the Scilab code into wxMaxima, currently I'm writing the tf2ss() function an inside that function the balanc() function is called, the problem is that I couldn't find the code for that function in Scilab installation directory, I've searched info in books and papers but every example starts with the Ab matrix given as an input to the problem, Scilab instead has the option to have as an input only the A matrix and it calculates the Ab and X matrices, so, I need help to make this function exactly as Scilab has it programmed to been able to compare all the steps that I'm doing.
Finally, wxMaxima has a function to calculate similarity transforms but it don't have the same output as Scilab what it means to me that they uses different criteria to calculate the similarity transform.
Note: I've tried to make the calculations in wxMaxima to have Ab and X matrices as elements with variables but the system of equations remains with too many variables and couldn't be solved.
Thanks in advance for the help in doing this.
In Scilab balanc() is hard-coded and based on LAPACK's dgebal (see the Fortran source at Netlib). In the algorithm the operations are quite simple (computing inf and 2-norms, swaping columns or rows of a matrix), maybe this could easily translated ?
A more readable version of the algorithm can be found on page 3 (Algorithm 2) of the following document: https://arxiv.org/abs/1401.5766.
Here is a Scilab implementation of Algorithm 3:
function [A,X]=bal(Ain)
A = Ain;
n = size(A,1);
X = ones(n,1);
β = 2; // multiply or divide by radix preserves precision
p = 2; // eventually change to 1-norm
converged = 0;
while converged == 0
converged = 1;
for i=1:n
c = norm(A(:,i),p);
r = norm(A(i,:),p);
s = c^p+r^p;
f = 1;
while c < r/β
c = c*β;
r = r/β;
f = f*β;
end
while c >= r*β
c = c/β;
r = r*β;
f = f/β;
end
if (c^p+r^p) < 0.95*s
converged = 0;
X(i) = f*X(i);
A(:,i) = f*A(:,i);
A(i,:) = A(i,:)/f;
end
end
end
X = diag(X);
endfunction
On this example the above implementation gives the same balanced matrix:
--> A=rand(5,5,"normal"); A(:,1)=A(:,1)*1024; A(2,:)=A(2,:)/1024
A =
897.30729 -1.6907865 -1.0217046 -0.9181476 -0.1464695
-0.5430253 -0.0011318 -0.0000356 -0.001277 -0.00038
-774.96457 3.1685332 0.1467254 -0.410953 -0.6165827
155.22118 0.1680727 -0.2262445 -0.3402948 1.6098294
1423.0797 -0.3302511 0.5909125 -1.2169245 -0.7546739
--> [Ab,X]=balanc(A)
Ab =
897.30729 -0.8453932 -32.694547 -14.690362 -9.3740507
-1.0860507 -0.0011318 -0.0022789 -0.0408643 -0.0486351
-24.217643 0.0495083 0.1467254 -0.2054765 -1.2331655
9.7013239 0.0052523 -0.452489 -0.3402948 6.4393174
22.23562 -0.0025801 0.2954562 -0.3042311 -0.7546739
X =
0.03125 0. 0. 0. 0.
0. 0.015625 0. 0. 0.
0. 0. 1. 0. 0.
0. 0. 0. 0.5 0.
0. 0. 0. 0. 2.
--> [Ab,X]=bal(A)
Ab =
897.30729 -0.8453932 -32.694547 -14.690362 -9.3740507
-1.0860507 -0.0011318 -0.0022789 -0.0408643 -0.0486351
-24.217643 0.0495083 0.1467254 -0.2054765 -1.2331655
9.7013239 0.0052523 -0.452489 -0.3402948 6.4393174
22.23562 -0.0025801 0.2954562 -0.3042311 -0.7546739
X =
1. 0. 0. 0. 0.
0. 0.5 0. 0. 0.
0. 0. 32. 0. 0.
0. 0. 0. 16. 0.
0. 0. 0. 0. 64.