Compute (skew)$\mu$ lower bounds for real perturbations.
Description
This routine computes (skew)$\mu$ lower bounds for an interconnection between a nominal LTI system $M(s)$ and a blockdiagonal operator $\Delta=diag\left(\Delta_1,\dots,\Delta_N\right)$, where $\Delta_i=\delta_iI_{n_i}$ and $\delta_i$ is real. At each point of a rough frequency grid, a power algorithm is run on a regularized problem and a series of LP problems is solved starting from the real part $\Delta_R$ of the resulting destabilizing perturbation (which usually brings the interconnection close to instability). The considered stability region can be bounded either by the imaginary axis or by a truncated sector.
Note: The aim is not to compute a (skew)$\mu$ lower bound as a function of frequency, but to compute a lower bound over the whole frequency range. The frequency at which the interconnection becomes unstable is left free during the optimization and can be very far from the initial one. Thus, a rough initial grid is usually sufficient to capture the most critical frequencies.
Syntax
[lbnd,wc,pert,tab]=mulb_nreal(sys,blk{,options})
Input arguments
The first two input arguments are mandatory:
sys  LTI object describing the nominal system $M(s)$. 
blk  Matrix defining the structure of the blockdiagonal operator $\Delta=diag\left(\Delta_1,...,\Delta_N\right)$. Its first 2 columns must be defined as follows for all $i=1,...,N$:
A third column can be used to specify skew uncertainties:
If none, 
The third input argument options
is an optional structured variable with fields:
freq  Frequency interval $\Omega$ in rad/s on which the lower bound is to be computed. The default value is options.freq=[0 10*max(abs(eig(sys)))] . 
grid  Frequency grid in rad/s ($1\times m$ array). It can also be a negative integer $m$, in which case a $m$point grid is automatically generated from $\Omega$. The default value is options.grid=10 . 
sector  Vector $[\alpha\ \xi]$ or $[\alpha\ \xi\ \omega_c]$ characterizing the considered truncated sector (see display_sector for a complete description). The default value is options.sector=[0 0] , which means that the left half plane is considered. 
target  The algorithm is interrupted if a lower bound is found which is greater than options.target . The default value is options.target=Inf . 
epsilon  A (skew)$\mu$ lower bound $\tilde{\mu}$ is acceptable if a perturbation $\tilde{\Delta}$ of norm 1/$\tilde{\mu}$ and a frequency $\omega_c$ are found, which satisfy $\mbox{det}(IM(j\omega_c)\tilde{\Delta})<\epsilon$, where $M(j\omega_c)$ is the frequency response of $M(s)$ at $j\omega_c$ (or at the corresponding point on the boundary of the truncated sector defined by options.sector ). The default value is options.epsilon=1e8 . 
trace  Trace of execution. The default value is options.trace=1 . 
warn  Warnings display. The default value is options.warn=1 . 
method  Type of LP problem:
The default value is

reg  Regularization to ensure convergence of the power algorithm. The default value is options.reg=0.05 , i.e. 5% of complex uncertainty is added. 
improve  If nonzero, additional computations are performed to find more accurate $\mu$ lower bounds. The default value is options.improve=0 . This optional parameter can only be used for classical $\mu$ problems. 
Output arguments
lbnd  Best (skew)$\mu$ lower bound on the considered frequency interval $\Omega$. 
wc  Frequency $\omega_c$ in rad/s for which lbnd has been computed. 
pert  Associated destabilizing perturbation $\tilde{\Delta}$: unless lbnd=0 , $\mbox{det}(IM(\omega_c)\tilde{\Delta})\approx 0$, where $M(\omega_c)$ is the frequency response of $M(s)$ at $j\omega_c$ (or at the corresponding point on the boundary of the truncated sector defined by options.sector ). 
tab  Structured variable with fields:

Example
System with 46 states and 20 nonrepeated real parametric uncertainties:load flexible_airplane
options.sector=[0.05 0.01];
[lbnd,wc,pert,tab]=mulb_nreal(sys,blk,options);lbnd
display_sector(options.sector,lft(sys,pert));
check_stab(lft(sys,pert),options.sector)
One of the eigenvalues of the interconnection lies just outside the considered truncated sector.abs(det(eye(20)calc_freq_resp(sys,wc,options.sector)*pert))
The value of $\mbox{det}(IM(\omega_c)\tilde{\Delta})$ is almost equal to 0.
See also
mulb
mulb_1real
mulb_mixed
hinflb_real
gen_grid
display_sector
Reference
[1]  G. Ferreres and JM. Biannic, "Reliable computation of the robustness margin for a flexible transport aircraft", Control Engineering Practice, vol. 9, pp. 12671278, 2001. 