The software RARL2 computes a stable rational L2-approximation of specified order n to a given stable p x m matrix function. We call stable matrix function a matrix whose entries belong to Hardy space H2 of the exterior of the disk. Namely, its entries are L2 functions on the unit circle that admit analytic continuation in the COMPLEMENT of the unit disk. This terminology comes from system theory: a stable matrix function is the transfer function of a multivariable discrete-time LTI system stable in the L2-Linfty sense.

The software RARL2 uses state-space representations of transfer functions which best suit numerical computations. A rational transfer function can be written in the form F(z)= D+C(z I-A)^{-1}B where A,B,C,D are matrices. The matrix A is square and is called the dynamics matrix. This (non unique) representation relates to the differential equations satisfied by the underlying dynamical system. The system is stable iff the eigenvalues of A belong to the unit disk. The 4-tuple (A,B,C,D) is called a realization of the system. If minimal, which means that the size of A is minimal among all admissible realizations, then the size of A is called the order or the McMillan degree of the system. The order of a scalar stable function p(z)/q(z) is precisely the polynomial degree of q(z) as expected. For more details on system and control theory, see [Kailath, Linear Systems, 1980]. The Control Toolbox Matlab uses both state-space models (called ss) and transfer functions models (called tf) and provides functions to convert one into the other.

The stable transfer function to be approximated can be given in one of the following forms
1) a realization (A,B,C,D)
2) its Fourier coefficients
3) pointwise values on the unit circle

Note that RARL2 requires the presence of the OPTIM and CONTROL toolboxes in order to work.

Content of the directory rarl2


boplib: the routines for a parametrization of stable allpass systems by means of balanced observable pairs

arl2lib: the routines for the use of “bop” parametrization for L2 rational approximation. It contains the functions arl2.m and rarl2.m which perform rational approximation from user data.

data: contain the matlab structure files (‘.mat’) of some selected examples listed below.


  • runcoeffHF : computes an approximation at order 8 of a 2×2 stable matrix function given by 501 matrix Fourier coefficients. These data comes from an application to the identification of microwave filters used in spacial telecommunications.
  • runfctKL08 : computes local approximants until order 9 of a scalar function of order 13. This example comes from [KL,2008].
  • runsampleKL08 : computes an approximation at order 8 of a the same function from 50 equidistant pointwise values on the circle.
  • runsyscontGlover : computes an approximation at order 8 of a 2×2 function of order 12. This example comes from [Glover,1984].
  • runsampleGlover : computes an approximation at order 8 of a the same function from 50 equidistant pointwise matrix values on the circle.
  • runsysMOHP02 : computes local approximants until order 3 of the function of order 1/z-1/z^3. At order 1, this real function possesses tree complex local approximants; only one of them is real and it is not the global minimum.

To handle your own data, you may use either the function arl2.m (simple search from an initial point which can either be provided by the user or automatically generated) or the function rarl2.m (iterative search starting from degree 1 and proceeding inductively; it requires no initial condition and is more likely to find the best approximant in case there are several local minima but is computationally heavier)

First add the required paths:


1) Enter the data structure.


data: matlab structure of type ‘coeff’; MIMO p x m system (the function to be approximated is p x m matrix valued, the Fourier coefficients are p x m complex matrices; the case of a scalar valued function requires setting p=m=1). The parameter N is the number of Fourier coefficients used to describe the fonction to be approximated, when the Fourier coefficient description is used (and not the realization description).

Latex formula

valInf(i,j)= F0(i,j)
userData.type = ‘coeff’
userdata.valinf = valInf (a pxm array)
userdata.value = CoeffPlus (a pxmxN array)

You may for example load the ‘coeff*.mat’ file available: ‘coeffHF’ :

load data/coeffHF

userData =
type: 'coeff'
value: [2x2x501 double]
valinf: [2x2 double]


data: matlab structure of type ‘sys’; MIMO p x m system.

F(z)=D+C(z I-A)^{-1}B

userData.type = ‘sys’
userData.value = a ss matlab structure
userData.value.a = A
userData.value.b = B
userData.value.c = C
userData.value.d = D

You may for example load one of the ‘.mat’ file available. For example load the system KL08 :

>> load data/sysKL08
>> userData

userData =

type: ‘sys’
value: [1×1 ss]

If your function is given as the quotient of two polynomials, then you may use the functions tf and ss (Control Toolbox) to convert the data into state space representation (see runfctKL08)

The following sequence produces a Nyquist plot of a state space model
>> [p,m] = size(theSys);
>> SSfig(p, m);
>> SSplot(theSys,’b’);
SSfig and SSplot are functions of the library arl2lib. Alternatively the Control Toolbox functions ‘nyquist’ and ‘bode’ can be used.


data: matlab structure of type ‘sample’; pointwise (p x m) matrix values on the unit circle: F(exp(i*theta(k))), k = 1,…,K.

You may for example load the ‘sample*.mat’ file available: ‘sampleKL08’ :

>> load data/sampleKL08
>> userData

userData =

type: ‘sample’
value: [1x1x50 double]
theta: [1×50 double]

2) Specify the rtype : ‘real’ or ‘complex’ approximation, according wether you have real Fourier coefficients and look for a best approximant with real Fourier coefficients as well, or if you want an approximant with possibly complex coefficients.

WARNING: the best complex rational approximant to a function with real Fourier coefficients needs not have real Fourier coefficients.

>> arl2_options=arl2Options(‘complex’)

to compute a best complex function approximating a given complex function.

>> arl2_options=arl2Options(‘real’)

if your function is real F: R^m –> R^p (i.e. possesses the conjugate symmetry) and you are looking for real approximants only.

WARNING: If your function does not possess conjugate symmetry, then real approximation does not work !

3) Specify optim_option, the options of the matlab function ‘fmincon’

Using the function optimset. If not specified, default values will be used. A convenient set of options is

>> optim_options = optimset(… …
‘MaxFunEval’, 10000, …
‘MaxIter’, 1000, …
‘Display’, ‘off’, …
‘TolFun’, 1.e-12, …
‘TolX’, 1.e-12, …
‘GradObj’, ‘on’ , …
‘DerivativeCheck’, ‘off’ …

The meaning of the above flags is of technical nature within the optimizing routine ‘fmincon’ and should be looked up in the documentation thereof. In particular, the termination tolerance on the criterion function value ‘TolFun’ should be carefully chosen. See the examples.

4) Choose your approximation method

arl2.m : simple approximation from a given initial point

Input arguments
userData.value: the matlab structure defined at step 1)
sys0: an initial point given by an ss matlab structure ‘sys0’
arl2_options: step 2)
optim_options: step 3)

The degree n of sys0 will be the order of the approximation.
An initial system can be randomly chosen using the RARL2 function ‘bopCreateStableSystem’
RARL2 also provides a function ‘fc2ss’ to compute an initial ‘system’ using Kung’s truncation method.

a) if userData are of ‘coeff’ type
>> n=8;
>> [sys0, Err] = fc2ss(userData.value,userData.valinf,n);

Err is the relative squared L2 error

b) if userData are of ‘sys’ type

>> [a,b,c,d]=ssdata(userData.value);
>> CoeffPlus=myImpulse(userData.value,50);
>> [sys0, Err] = fc2ss(CoeffPlus,d,n);

c) if userData are of ‘sample’ type

>> [p,m,NC] = size(userData.value);
>> n=8 –> approximant degree
>> valinf = zeros(p,m);
value = zeros(p,m,NC);
>> for lig=1:p
for col=1:m
f = userData.value(lig,col,:);
t = userData.theta;
sample_fc = struct(‘t’,t,’f’,f);
fc = sample2fc(sample_fc, 1, NC,’rtype’); –> rtype=’real’ or ‘complex’
>> [sys0, Err] = fc2ss(value,valinf,n);

Then GO

>> report=arl2(userData,sys0,arl2_options,optim_options)

report is a matlab structure which contains some of the output arguments of fmincon and the result of the approximation

report =

iterations: 27
funcCount: 500
lssteplength: 1
stepsize: 3.3392e-07
algorithm: ‘medium-scale: SQP, Quasi-Newton, line-search’
firstorderopt: 7.5971e-08
constrviolation: -0.9998
message: [1×776 char]
finalValue: 2.0129e-06 –> square of L2 error norm
options: [1×1 struct]
CSP: [1×1 struct]
execTime: 11.7900
sys: [2×2 ss] –> the approximant (ss structure)


RARL2 recursively runs arl2: when processing degree k+1, a number of maxIP x p initial guesses
are built from the best k-degreed order approximant to arl2.

Input arguments
userData: the matlab structure defined at step 1)
n: required approximant degree
maxIP: default is 2 in the case ‘real’ and 4 in the case ‘complex’
arl2_options: step 2)
optim_options: step 3)

>> report=RARL2(userData,n,maxIP,arl2_options,optim_options)

report is a matlab structure which contains all the approximants from degree 1 to target degree n.

data: [1×1 struct]
n: 3
maxIP: 4
exectime: 44.5800
LM0: {[1×1 struct]}
LM: {{1×3 cell} {1×3 cell} {1×1 cell}}
J0: 2.0000
message: ‘data & approximants’
sys: [1×1 ss]
finalValue: 2.7698e-15

The data and local minima at each degree can be plotted using the RARL2 function ShowReport. These plots are Nyquist diagrams, which allows for a quick estimation of the quality of the approximation.

>> ShowReport(report)

[KL,2008]: A. Lindquist and J. Karlsson, Stability-preserving rational approximation subject to interpolation constraints, IEEE Trans. Automatic Control 53(7) , 1724–1730 (2008).
[Glover, 1984]: K. Glover, All optimal Hankel-norm approximations of linear multivariable systems and their error bounds, Int. J. Control 39, 1115-1193 (1984).
[MOHP, 2002]: J.P. Marmorat, M. Olivi, B. Hanzon, R. Peeters, Matrix rational H2-approximation: a state-space approach using Schur parameters , CDC02 (Las Vegas)

Copyright©2002 INRIA-Mines ParisTech
Jean-Paul MARMORAT & Martine OLIVI

email: Martine.Olivi@sophia.inria.fr

Web http://www-sop.inria.fr/teams/apics/RARL2/rarl2.html

Comments are closed.