h2invsyslin — Stable inverse with minimal weighted H2 norm.
SYSG = h2invsyslin(SYSH [,opts])
[SYSG,SYSA,SYSB] = h2invsyslin(SYSH [,opts])
Discrete-time state-space realization of H(z) (a syslin list). The matrix SYSH.D should be of full rank. The initial state SYSH.x0 may be non-zero.
Discrete-time state-space realization of G(z) (a syslin list).
Discrete-time state-space realization of A(z) (a syslin list).
Discrete-time state-space realization of B(z) (a syslin list).
A scilab struct which may contain additional options. Optional parameter, default value is struct().
Delay in the output of the inverse G(z) (a natural). Optional parameter, default value is zero.
Discrete-time state-space realization of W(z) (a syslin list), which must have as many outputs as H(z). SYSW has to be stable. Optional parameter, default value is syslin("d",[],[],[],eye(,)).
The routine computes a stable state-space system system SYSG such
that its transfer function G(z) satisfies G(z)*H(z)=I/z^L
and ||G*W||=min
, where L is short for opts.delay. The
complete set of minimum norm inverses of H(z) is given by
A(z)+K(z)*B(z)
, where K(z) is an arbitrary stable and
causal transfer function.
Assumptions: H(z) should have full normal rank. The finite zeros of SYSH, i.e., the complex z where the rank of the matrix pencil [SYSH.A-z*eye() SYSH.B;SYSH.C SYSH.D] drops below its normal rank, must be located inside the open unit disc |z|<1. There also should be no zero at infinity, i.e., SYSH.D should be of full rank. If weighting is used, i.e., a custom opts.weight is specified, then [I-H(z)*pinv(H(z))]*W(z) should have constant rank on the unit circle.
Note that if H(z) is wide and not tall, optimal right-inverses are
computed instead of left-inverses, i.e., H(z)*G(z)=I/z^L
and
||W*G||=min
. The set of optimal right-inverses becomes
A(z)+B(z)*K(z)
. The last assumption changes to
W(z)*[I-pinv(H(z))*H(z)] having constant rank on the unit circle.
// Example 1: Minimal norm left inverse of the // transfer function [1/z;1-1/z] with delay three z = poly(0,"z"); // define symbolic "z" variable H = [1/z ; 1-1/z]; // define transfer function SYSH = tf2ss(H); // convert H to state-space SYSH.dt = "d"; // mark as discrete-time opts = struct("delay",3); // define the delay SYSG = h2invsyslin(SYSH,opts); // compute left inverse G = ss2tf(SYSG); // convert inverse into rational // matrix disp("Left Inverse G of H"); disp(clean(G)); // show results disp("H2 Norm of G:"); disp(dh2norm(SYSG)); disp("G*H:"); disp(clean(G*H)); // Example 2: Recovery of input from output // for a random system p = 2; q = 4; N = 5; n = N*p; // create random 5 tap FIR A = [zeros(p,n) ; eye(n-p,n-p) zeros(n-p,p)]; // system with 2 inputs and B = [eye(p,p) ; zeros(n-p,p)]; // 4 outputs and random initial C = rand(q,n); // state D = rand(q,p); x0 = rand(n,1); SYSH = syslin("d",A,B,C,D,x0); K = 100; // define no. of input vectors U = ones(p,K); // create input vectors L = 10; // L is the delay Y = dsimul(SYSH,[U zeros(p,L)]); // apply SYSH to the input opts = struct("delay",L); SYSG = h2invsyslin(SYSH,opts); // compute left inverse U2 = dsimul(SYSG,Y); // apply SYSG to output of SYSH disp("Residual:"); disp(norm(U-U2(:,L+1:$))); // check that U2 is the delayed // input, i.e., that the // residual is approx. zero // Example 3: Inverse with minimal weighted H2 norm p = 2; q = 3; N = 5; n = N*p; // create random 5 tap FIR A = [zeros(p,n) ; eye(n-p,n-p) zeros(n-p,p)]; // system with 2 inputs and B = [eye(p,p) ; zeros(n-p,p)]; // 3 outputs and random initial C = rand(q,n); // state D = rand(q,p); SYSH = syslin("d",A,B,C,D); p = 5; q = 3; N = 5; n = N*p; // create random 5 tap FIR A = [zeros(p,n) ; eye(n-p,n-p) zeros(n-p,p)]; // weight with 5 inputs and B = [eye(p,p) ; zeros(n-p,p)]; // 3 outputs and random initial C = rand(q,n); // state D = rand(q,p); SYSW = syslin("d",A,B,C,D); opts = struct("delay",3); // define the delay SYSG = h2invsyslin(SYSH,opts); // compute optimal inverse // without weight disp("G*H:"); disp(clean(ss2tf(SYSG*SYSH))); // show results disp("H2 norm of G*W:"); disp(dh2norm(SYSG*SYSW)); opts = struct("delay",3,"weight",SYSW); SYSG = h2invsyslin(SYSH,opts); // compute inverse with weight disp("G*H:"); disp(clean(ss2tf(SYSG*SYSH))); // show results // weighted norm should be // smaller than before! disp("H2 norm of G*W:"); disp(dh2norm(SYSG*SYSW));