System L1-norm



01 May 2013 (Updated )

Calculate L1-norm of impulse response of continuous-time SISO LTI system

%   function [L1norm,err,U,L,n,tol]=l1norm(G,tol);
%	l1norm(G)
% 	Calculate L1-norm (Rutland & Lane's algorithm)
% 	of the impulse response of a proper SISO LTI system with
%	minimal state space realization G=(A,B,C,0)
%	[norm,err,U,L,n,tol]=l1norm(G)
%	Additional outputs:
%	err - error estimate
%	U - upper bound on norm
%	L - lower bound on norm
%	niter - number of iterations
%	tol = -alpha * 1e-3 where alpha is maximum
%	      real part of the eigenvalues of A
%	l1norm(G,tol)
% 	Calculate L1-norm by Rutland & Lane's algorithm
% 	of the impulse response of a SISO LTI system with
%	minimal state space realization (A,B,C) to a
%	tolerance tol
%   (default tol=-alpha/1000 where alpha is the abscissa of stability)
%	l1norm(G,tol,maxiter)
% 	Calculate L1-norm by Rutland & Lane's algorithm
% 	of the impulse response of a SISO LTI system with
%	minimal state space realization (A,B,C) to a
%	tolerance tol with a maximum number of of iterations niter 
%   (default maxiter=20) 
%   NOTE that memory problems can occur for niter>20, so increase with
%        caution
%	See Rutland N.K. & Lane P.G, "Computing the 1-norm of the impulse 
%        response of linear time-invariant systems", 
%        Systems and Control Letters, Volume 26, Number 3, pp. 211-221 

% (c) JF Whidborne 28/4/95 
%     modified jfw 1/5/2013, 6/10/14, 13/11/14

function [L1norm,err,U,L,tol,niter]=l1norm(G,tol,maxiter)

Gss=ss(G); % ensure state space realization

% check SISO
[ny,nu] = size(Gss.d);
if ny~=1||nu~=1
    error('Error: System not SISO')

% Check strict properness
if Gss.d~=0
    error('Error: System not strictly proper')

A=Gss.a;B=Gss.b;C=Gss.c; % get state matrices/vectors
% abscissa of stability alpha
if alpha>0, error('unstable system');end
if cond(A)>1e6, % 13/11/14
    warning('ill conditioned system - results may be inaccurate - use balreal to obtain a balanced realization');

if nargin<3, maxiter=20; end% maximum iterations
if nargin<2, tol=-alpha*.001; end% tolerance
[ns,~] = size(A); % state dimension
x0 = zeros(1,ns);% initial condition
G0=-C*inv(A)*B;% final steady state value (step response)

% 1st step T=0, N=0

[~,nu_u_Nplus1] = fminbnd(@normtail,0,-alpha,...
%fplot(@(sig) normtail(sig,A,B',C,ns),[0,-alpha]);

nu_u(1)= nu_u_Nplus1;%  weighted H2 norm
U=nu_u(1); % upper bound
nu_l(1)=abs(G0); % steady state for step response
L=nu_l(1); % lower bound

niter=1;N=1;T=-5/alpha; T_flag=1; % estimate of required time T

while((((U-L)/L/2)>=tol) && niter<maxiter)% increase maxiter for accuracy if memory available !!
    niter=niter+1; % iteration number
    h=T/N; % step length
    [Phi, Gamma] = c2d(A,B,h); % get difference equation
    % evaluate nu_l(n) lower bound
    y = ltitr(Phi,Gamma,[1;zeros(N,1)],x0)*C.';% get impulse response
    nu_l_k=abs(y(2:N+1));% lower bounds in intervals
    if T_flag, % time has changed - get nu_l_{N+1} error
        %y1 = ltitr(Phi,Gamma,ones(N+1,1),x0)*C.';% get 1 response
        %nu_l_Nplus1= abs(G0- y1(N+1));
        nu_l_Nplus1= abs(G0- sum(y));%sum(y) is step response at t=T
    end% if
    sum_nu_l_k=sum(nu_l_k); % lower bound sum <=T
    nu_l(niter) = sum_nu_l_k + nu_l_Nplus1;% eqn(8)
    % evaluate nu_u(n) upper bound
    x = ltitr(Phi,Gamma,zeros(N+1,1),B);% state impulse response, get value of last point
    CTC=C'*C; W=lyap(A',CTC-Phi'*CTC*Phi);% Lyapanov eqn(36)
    %nu_u_k=(h*sum(((x(1:N,:)*W).*x(1:N,:))')).^.5;% unrefined upper bounds
    %sum_nu_u_k=sum(nu_u_k); % unrefined upper bound sum <=T
    % now evaluate refined upper bounds
    int_diff=sqrt(h*sum(((x(1:N,:)*AWA).*x(1:N,:))'));% eqn (41)
    n_rub=find((max(e_abs(1:N),e_abs(2:N+1)))<=int_diff);% intervals which do not have exact lower bound
    nu_u_k(n_rub)=sqrt(h*sum(((x(n_rub,:)*W).*x(n_rub,:))'));% refined upper bounds
    if T_flag, % time has changed - get nu_u_{N+1} error
%        [sig_min,OPTIONS]=fmin('normtail',0,-alpha,OPTIONS1,A,x(N+1,:),C,ns);
        [~, nu_u_Nplus1] = fminbnd(@normtail,0,-alpha,...
        T_flag=0; % reset flag
    end% if
    sum_nu_u_k=sum(nu_u_k); % upper bound sum <=T
    nu_u(niter) = sum_nu_u_k + nu_u_Nplus1;% eqn(11)
    % calculate lower bound L
    % calculate upper bound U
    %double T criterion - change final time and set flag
err=(U-L)/L/2;% error estimate
L1norm=(U+L)/2; % mean of upper & lower bounds

function nu=normtail(sig,A,x,C,ns)
% minimising function used by l1norm.m - eqn(39)
% 1/2\sigma \int_0^\infty |\exp(\sigma t) e(t+T,\delta)|^2 dt

nu= ((x*W*x')/2/sig)^(0.5);%need square root

Contact us