%% Análisis de márgenes de estabilidad con modelos inciertos que tienen rango asimétrico
% © 2020 Antonio Sala Piqueras, Universitat Politècnica de València. Todos los derechos reservados.
% Presentación en vídeo: http://personales.upv.es/asala/YT/V/asimrc.html
% Este código funcionó correctamente en Matlab R2020a
% Objetivo: Las incertidumbres asimétricas son "transformadas" a la imagen de una bola simétrica de tamaño 1 (ver "Decomposing Uncertain Objects" en la documentación de Matlab) que es la necesaria en la teoría. Este código proporciona un ejemplo para comprender dicha transformación.
%% Modelo con parámetros inciertos

p=ureal("P",.5,"PlusMinus",[-.25,+.1]);
q=ureal("Q",5,"PlusMinus",[-1,+3]);
s=tf('s');
G=1/(s^2+q*(s+p)+5);

%% Análisis de resultados de robstab
% El robust control toolbox hace una transformación de un parámetro \rho de la forma F(\rho)=(a+b\rho)/(c+d\rho). Los coeficientes a, b, c y d son calculados de modo que el valor nominal se mapea a 0, y los valores extremos del rango a -1 y +1. Veamos el ejemplo:

opts=robOptions('Display','on','Sensitivity','on');
[rs,wcu]=robstab(G,opts)
pole(usubs(G,wcu))

%% Interpretación gráfica de los márgenes

dibujamargenes(p,[-.35 .75],rs,wcu.P)
dibujamargenes(q,[3 21],rs,wcu.Q)

%% Conclusiones
% Los resultados teóricos todos hacen referencia a "bolas" de incertidumbre de radio unidad.
% Los rangos [LB,N,UB] no simétricos son transformados al rango [-1,0,1], y los márgenes de estabilidad robstab deben ser interpretados en ese rango, que puede significar otro rango diferente en unidades físicas (mucho más asimétrico).
% Los comandos de Matlab encargados de la transformación y su inversa son actual2normalized, normalized2actual.
%
%

function dibujamargenes(q,rango,rs,peorcaso)
  AV = linspace(rango(1),rango(2),250); %AV : rango en unidades físicas para plot
  NV = actual2normalized(q,AV); %NV: unidades normalizadas correspondientes
  plot(NV,AV,0,q.NominalValue,'ro',-1,q.Range(1),'bo',1,q.Range(2),'bo','Linewidth',2) 
  grid, xlabel('Valor Normalizado'), ylabel('Valor físico de q')
  xline(-1,'-.g','Lím.inf');xline(1,'-.g','Lím.sup'); %axis tight
  xline(-rs.LowerBound,'-.r','robstab');xline(rs.LowerBound,'-.r','robstab'); 
  xline(0,'-.b','Nominal');
  yline(peorcaso,'--k','Peor Caso')
end

%