mu, muunwrap, randel, unwrapd, unwrapp (Mu Analysis and Synthesis Toolbox)

Mu Analysis and Synthesis Toolbox

mu, muunwrap, randel, unwrapd, unwrapp

Compute upper and lower bounds for the complex and mixed (real and complex) structured singular value (referred to as mixed µ) of a VARYING/CONSTANT matrix

Syntax

[bnds,dvec,sens,pvec,gvec] = mu(matin,blk,options);
[dl,dr,gl,gm,gr,pert] = muunwrap(dvec,gvec,pvec,blk);
[dl,dr,gl,gm,gr] = muunwrap(dvec,gvec,blk);
[dl,dr] = muunwrap(dvec,blk);
pert = randel(blk,nrm,opt);
[dl,dr] = unwrapd(dvec,blk);
pert = unwrapp(pvec,blk);

The functions associated with mixed µ are:


`mu`	general complex and mixed computation
`muunwrap`	extract block-diagonal D- and G- scalings from row vectorscontaining the scalings
`randel`	create random block structured matrix
`unwrapd`	extract block-diagonal D-scalings from row vectorcontaining the scalings
`unwrapp`	extract block-diagonal perturbation from row vector containing the perturbation

Description
Input arguments:

matin
A CONSTANT or VARYING matrix

blk
An array that describes the perturbation block structure. Its size is nblk x 2, where nblk is the total number of blocks in the perturbation structure. The ith row of blk defines the dimensions of the ith perturbation block. If blk(i,:) = [-r 0], then the ith block is an r x r repeated, diagonal real scalar perturbation, while if blk(i,:) = [r 0], then the ith block is an r x r repeated, diagonal complex scalar perturbation, and if is blk(i,:) = [r c], then the ith block is an r x c complex full-block perturbation. If blk is omitted, the default is all 1 x 1 complex blocks, and results in an error if matin is not square.

options

An optional character string describing the desired computations. It can consist of the following characters:
'l' - compute lower bound using a power iteration
't' - use more iterations in the lower bound
'R' - start power iteration with RANDOM vectors
'Rj' - restart lower bound j times with RANDOM vectors where j is an integer between 1 and 9
'u' - compute upper bound using a balanced/LMI technique
'c' - compute upper bound to greater accuracy'C'
'C' - compute tighest upper bound (may be slow)
'f' - compute a fast but crude upper bound
'r' - restart computation at EACH independent variable
's' - suppress progress information
'w' - suppress warnings
'L' - compute only the lower bound
'U' - compute only the upper bound

`matin`	A CONSTANT or VARYING matrix
`blk`	An array that describes the perturbation block structure. Its size is `nblk` x 2, where `nblk` is the total number of blocks in the perturbation structure. The ith row of `blk` defines the dimensions of the ith perturbation block. If `blk(i,:) = [-r 0]`, then the ith block is an r x r repeated, diagonal real scalar perturbation, while if `blk(i,:) = [r 0]`, then the ith block is an r x r repeated, diagonal complex scalar perturbation, and if is `blk(i,:) = [r c]`, then the ith block is an r x c complex full-block perturbation. If `blk` is omitted, the default is all 1 x 1 complex blocks, and results in an error if `matin` is not square.
`options`	An optional character string describing the desired computations. It can consist of the following characters: `'l'` - compute lower bound using a power iteration `'t'` - use more iterations in the lower bound `'R'` - start power iteration with RANDOM vectors `'Rj'` - restart lower bound `j` times with RANDOM vectors where `j` is an integer between 1 and 9 `'u'` - compute upper bound using a balanced/LMI technique `'c'` - compute upper bound to greater accuracy`'C'` `'C' -` compute tighest upper bound (may be slow) `'f'` - compute a fast but crude upper bound `'r'` - restart computation at EACH independent variable `'s'` - suppress progress information `'w'` - suppress warnings `'L'` - compute only the lower bound `'U'` - compute only the upper bound

The default value of options is 'lu', meaning that a lower bound will be computed using the power method, Young and Doyle 1990 and Packard et al. 1988, and an upper bound will be computed, using the balanced/AMI technique, Young et al., 1992, for computing the upper bound from Fan et al., 1991.

Output arguments:

bnds
A 1 x 2 vector. If matin is VARYING, so is bnds, whereas if matin is a CONSTANT matrix, then bnds is CONSTANT. The first column of bnds contains an upper bound to mixed µ of matin, and the second column contains a lower bound to mixed µ.

dvec and gvec
Row vectors which contain the D and G scaling matrices that have produced the upper bound in bnds. dvec and gvec are the same data type as bnds and are stored as vectors to save memory. They can be unwrapped into the appropriate D and G matrices by using the command, muunwrap:

`bnds`	A 1 x 2 vector. If `matin` is VARYING, so is `bnds`, whereas if `matin` is a CONSTANT matrix, then `bnds` is CONSTANT. The first column of `bnds` contains an upper bound to mixed µ of `matin`, and the second column contains a lower bound to mixed µ.
`dvec and gvec`	Row vectors which contain the D and G scaling matrices that have produced the upper bound in `bnds`. `dvec` and `gvec` are the same data type as `bnds` and are stored as vectors to save memory. They can be unwrapped into the appropriate D and G matrices by using the command, `muunwrap`:

    [dl,dr,gl,gm,gr] = muunwrap(dvec,gvec,blk); 

The upper bound in bnds for a matrix M is a number  > 0 such that there are scaling matrices Dl, Dr, Gl, Gm, Gr (see Young et al., 1992, for details) satisfying



sens is a row vector which contains the sensitivity of

with respect to the values in Dl (and Dr). It is calculated in an ad-hoc manner, and is mainly used when fitting frequency varying Ds with rational functions via the routines in dkit, dkitgui, and autodkit.

pvec is a row vector containing a perturbation matrix that has the structure defined by blk. As with dvec and gvec, pvec is the same data type as matin and is stored as a vector. It can be unwrapped into the actual perturbation with the command, unwrapp:

pert = unwrapp(pvec,blk);

After being unwrapped, the perturbation matrix pert satisfies three conditions:

It has the block structure defined by blk;
The maximum singular value of pert is equal to the reciprocal of the lower bound in bnds (when the lower bound is not zero);
The matrix mmult(matin,pert) has an eigenvalue equal to 1 at each independent variable.

Examples
Suppose sys is a system matrix with four inputs and four outputs, and that it is stable. sys_g is a frequency response of sys.

is 4 1 x 1 perturbation blocks

```
blk = ones(4,2); 
```

% Calculate µ on frequency response

[bnds,dvec,sens,pvec] = mu(sys_g,blk);

% Unwrap the D scaling matrices

```
[dl,dr] = unwrapd(dvec,blk); 
```

% Generate scaled matrix

```
dmdi = mmult(dl,sys_g,minv(dr)); 
```

% Verify the upper bound

vplot('liv,m',vnorm(dmdi),sel(bnds,1,1))

% Unwrap the perturbation

```
actpert = unwrapp(pvec,blk); 
```

% Check that perturbation is correct structure

```
see(xtracti(actpert,1)) 
```

% Check lower bound and perturbation

vplot('liv,lm',sel(bnds,1,2),vnorm(actpert))

% Form M

```
mdel = mmult(sys_g,actpert); 
```

% Verify M has an eigenvalue at 1

```
veig(mdel)
```

Looking at the same frequency response sys_g with a mixed real/complex block structure:

is 4 1 x 1 perturbation blocks

% with the first 2 real, and the last 2 complex

```
blk = [-1 0; -1 0; 1 0; 1 0]; 
```

% Calculate mixed µ on frequency response

[bnds,dvec,sens,pvec,gvec] = mu(sys_g,blk);

% Unwrap the D and G scaling matrices

[dl,dr,gl,gm,gr] = muunwrap(dvec,gvec,blk);

% Generate scaled matrix

dmd = mmult(dl,sys_g,minv(dr)); 
oobdmd = veval('*',dmd,minv(sel(bnds,1,1))); 
oobdmdjg = msub(oobdmd,mscl(gm,j)); 
scall = madd(eye(4),mmult(gl,gl)); 
scalr = madd(eye(4),mmult(gr,gr)); 
scall = veval(,scall,-0.25); 
scalr = veval(",scalr,-0.25); 
scaledmat = mmult(scall,oobdmdjg,scalr);

% Verify the upper bound (scaledmat should have norm 1)

```
vplot('liv,m',vnorm(scaledmat)) 
```

% Unwrap the perturbation

```
pert = unwrapp(pvec,blk); 
```

% Check that perturbation is correct structure

```
see(xtracti(pert,1)) 
```

% Check lower bound and perturbation

vplot('liv,m',sel(bnds,1,2),minv(vnorm(pert)))

% Form M

```
mdel = mmult(sys_g,pert); 
```

% Verify M has an eigenvalue at 1

```
see(veig(mdel))
```

For more examples of computing bounds for µ, please refer to the "Computational Exercise with the mu Command" and "Computational Exercise with the mu Command -- Mixed Perturbations" sections in Chapter 4 as well as the robust multivariable control design examples in Chapters 6 and 7.

Algorithm
Peter Young and Matt Newlin helped write the mu program and supporting routines.

The lower-bound power algorithm is from Young and Doyle, 1990, and Packard et al. 1988.

The upper-bound is an implementation of the bound from Fan et al., 1991, and is described in detail in Young et al., 1992. In the upper bound computation, the matrix is first balanced using either a variation of Osborne's method (Osborne, 1960) generalized to handle repeated scalar and full blocks, or a Perron approach. This generates the standard upper bound for the associated complex µ problem. The Perron eigenvector method is based on an idea of Safonov, (Safonov, 1982). It gives the exact computation of µ for positive matrices with scalar blocks, but is comparable to Osborne on general matrices. Both the Perron and Osborne methods have been modified to handle repeated scalar and full blocks. Perron is faster for small matrices but has a growth rate of n³, compared with less than n² for Osborne. This is partly due to the MATLAB implementation, which greatly favors Perron. The default is to use Perron for simple block structures and Osborne for more complicated block structures. A sequence of improvements to the upper bound is then made based on various equivalent forms of the upper bound. A number of descent techniques are used which exploit the structure of the problem, concluding with general purpose AMI optimization (Boyd et al.), 1993, to obtain the final answer.

Reference
S. Boyd, and L. El Ghaoui, "Methods of centers for minimizing generalized eigenvalues," Linear Algebra and Its Applications, vol. 188-189, pp. 63-111, 1993.

M. Fan, A. Tits, and J. Doyle, "Robustness in the presence of mixed parametric uncertainty and unmodeled dynamics," IEEE Transactions on Automatic Control, vol. AC-36, pp. 25-38, 1991.

E. Osborne, "On preconditioning of matrices," Journal of Associated Computer Machines, vol. 7, pp. 338-345, 1960.

A. Packard, M. Fan and J. Doyle, "A power method for the structured singular value," Proc. of 1988 IEEE Conference on Control and Decision, pp. 2132-2137, December 1988.

M. Safonov, "Stability margins for diagonally perturbed multivariable feedback systems," IEEE Proc., vol. 129, Part D, pp. 251-256, 1982.

P. Young and J. Doyle, "Computation of with real and complex uncertainties," Proceedings of the 29th IEEE Conference on Decision and Control, pp. 1230-1235, 1990.

P. Young, M. Newlin, and J. Doyle, "Practical computation of the mixed problem," Proceedings of the American Control Conference, pp. 2190-2194, June, 1992.

See Also
blknorm, dypert, genmu, norm, vnorm, vrho, vsvd

mscl, sclin, sclout musynfit, musynflp, muftbtch