%Figure 8.15 (Table 8.1, mu-analysis)
clear all;close all;
G0 = [87.8 -86.4; 108.2 -109.6];
dyn = nd2sys(1,[75 1]); Dyn=daug(dyn,dyn); 
G=mmult(Dyn,G0);
% Inverse-based control
dynk=nd2sys([75 1],[1 1.e-5],0.7); Dynk=daug(dynk,dynk); 
Kinv=mmult(Dynk,minv(G0));
% Weights
wp=nd2sys([10 1],[10 1.e-5],0.5); Wp=daug(wp,wp); 
wi=nd2sys([1 0.2],[0.5 1]); Wi=daug(wi,wi);
% Generalized plant P
systemnames = 'G Wp Wi';
inputvar = '[ydel(2); w(2) ; u(2)]'; outputvar = '[Wi; Wp; -G-w]';
input_to_G = '[u+ydel]'; input_to_Wp = '[G+w]'; input_to_Wi = '[u]';
sysoutname = 'P'; cleanupsysic = 'yes'; sysic;
N = starp(P,Kinv);
omega = logspace(-3,3,61); Nf = frsp(N,omega);
% mu for RP
blk = [1 1; 1 1; 2 2]; [mubnds,rowd,sens,rowp,rowg] = mu(Nf,blk,'c'); 
muRP = sel(mubnds,':',1); pkvnorm(muRP)           % (ans = 5.7726)
% Worst-case weighted sensitivity (could also use skewmu-command)
[delworst,muslow,musup] = wcperf(Nf,blk,1);
mus=vunpck(vinterp(musup,[1 1;1 Inf],1)) % (musup = 44.92 for delta=1)
delworst=delworst
% mu for RS
Nrs=sel(Nf,1:2,1:2); [mubnds,rowd,sens,rowp,rowg]=mu(Nrs,[1 1; 1 1],'c');
muRS = sel(mubnds,':',1); pkvnorm(muRS)           % (ans = 0.5242)
% mu for NP (= max. singular value of Nnp)
Nnp=sel(Nf,3:4,3:4); [mubnds,rowd,sens,rowp,rowg]=mu(Nnp,[2 2],'c');
muNP = sel(mubnds,':',1); pkvnorm(muNP)           % (ans = 0.5000)
vplot('liv,m',muRP,muRS,'--',muNP,'-.');  %Figure 8.15
xlabel('Frequency');ylabel('ssv');
text(0.01,0.5,'RS');
text(20,0.8,'NP');
text(0.3,4.5,'RP');
