<< dh2norm Linear System Inversion Toolbox h8invsyslin >>

Linear System Inversion Toolbox >> Linear System Inversion Toolbox > h2invsyslin

h2invsyslin

Stable inverse with minimal weighted H2 norm.

Calling Sequence

SYSG = h2invsyslin(SYSH [,opts])
[SYSG,SYSA,SYSB] = h2invsyslin(SYSH [,opts])

Parameters

SYSH

Discrete-time state-space realization of H(z) (a list). The matrix SYSH.D should be of full rank. The initial state SYSH.x0 may be non-zero.

SYSG

Discrete-time state-space realization of G(z) (a list).

SYSA

Discrete-time state-space realization of A(z) (a list).

SYSB

Discrete-time state-space realization of B(z) (a list).

opts

A scilab which may contain additional options. Optional parameter, default value is struct().

opts.delay

Delay in the output of the inverse G(z) (a natural). Optional parameter, default value is zero.

opts.weight

Discrete-time state-space realization of W(z) (a list), which must have as many outputs as H(z). SYSW has to be stable. Optional parameter, default value is syslin("d",[],[],[],eye(,)).

Description

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.

Examples

// 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));

<< dh2norm Linear System Inversion Toolbox h8invsyslin >>