%% Pasividad + LFT: criterio del círculo
% © 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/pascricir.html
% Este código funcionó con Matlab R2020a
% Objetivo: Dada la interconexión
%
% evaluar cómo debe ser la respuesta en frecuencia de G para que la interconexión sea estable, para caso serie+paralelo de F.
% Tabla de Contenidos
% Revisión teórica
% Impedancias serie o paralelo
%     Caso 1:  exceso de pasividad input
%     Caso 2, exceso de pasividad output
% Caso 3, exceso de pasividad input+output (simultáneo)
%     Estabilidad robusta
%         Criterio del círculo (realim +)
% Conclusiones
% Funciones auxiliares
%% Revisión teórica
%
% H=lft(F,G)
% El Teorema de Pasividad dice que si \Delta es estrictamente pasivo y -H también (o estrict. output passive y observable), entonces la interconexión es estable. (Nyquist de H en lado izquierdo del plano complejo, parte real negativa)
%
% Mirando únicamente la parte de abajo, H=lft(F,G), m=Hn. Entonces n=H^{-1}m, con lo que se debe cumplir u=lft(H^{-1},F)\cdot y, esto es G=\left(lft\left(H^{-1},F\right)\right)^{-1} .
% Si Re(H)< 0, entonces Re(H^{-1})< 0 -- en efecto, si H=a+bj entonces H^{-1}=\frac{a}{a^2+b^2}-\frac{b}{a^2+b^2}j.
% Si F no depende de la frecuencia (el caso Fdinámico es más complicado), el lugar geométrico de la respuesta en frecuencia de G es:
% \bigcup_{\sigma>0,\gamma\in\mathbb R}\, \frac{1}{lft(-\sigma+\gamma j,F)}
%% Impedancias serie o paralelo
%% Caso 1:  exceso de pasividad input
% Ejemplo eléctrico: resistencia en serie con circuito pasivo.
%
% Las mismas ecuaciones tiene el sistema mecánico [entrada v, salida F]:
%  F=F_1+F_p, F_1=Rv
%
% m=y, u=n+\alpha y
% Si \int mn \,dt\geq 0, entonces \int uy\,dt=\int (n+\alpha y)y\,dt=\int (nm+\alpha y^2)\,dt\geq \alpha \int y^2\,dt

alpha=3;
F_serie=[0 1;1 alpha];
dibujatodo(F_serie)
hold on, plot(1/alpha,0,'rx'), hold off, axis([-1.8 1.8 -1.5 1.5])

%
%
%% Caso 2, exceso de pasividad output
% Ejemplo eléctrico: resistencia en paralelo con circuito pasivo.
%
% Las mismas ecuaciones tiene el sistema mecánico [entrada v, salida F]:
%
% v=v_1+v_p, v_1=F/R
% En general, consideramos:
% m=-\beta n+y, u=n
% Si \int mn \,dt\geq 0, entonces \int uy\,dt=\int n(m+\beta n)\,dt=\int (nm+\beta u^2)\,dt\geq \beta \int u^2\,dt

beta=1/4;
F_paralelo=[-beta 1;1 0];
figure
dibujatodo(F_paralelo)
hold on, plot(beta,0,'rx'), hold off
axis([-2 0.5 -2 2])

% G debe ser estable.
%%
%
%% Caso 3, exceso de pasividad input+output (simultáneo)
% Consideremos la conexión:
%
% La impedancia resistiva (parte real) mínima es R_{\min}=R_1 cuando z\to 0, y la máxima es R_{\max}=R_1+R_2 cuando z\to\infty.
% El componente desconocido z puede ser lineal o no lineal, estático o dinámico, mientras sea pasivo.
% Análogo mecánico: cambiando tensión a fuerza, intensidad a velocidad, tenemos F=F_p+F_1, F_1=R_1v, F_p=R_2v_2, v_2+v_p=v, siendo las velocidades entendidas como diferencia de velocidad entre extremos de los amortiguadores o del sistema mecánico pasivo desconocido. [entrada v, salida F]
%

Rmin=alpha, Rmax=1/beta+alpha

% Definiendo:                \beta=\frac{1}{R_2},         \alpha=R_1
% El cambio de variable (LFT) será:
% m=-\beta n+y, u=n+\alpha y
% \left(\begin{array}{c}m\\u \end{array}\right)=\underbrace{\left(\begin{array}{cc}-\beta &1\\1 &+\alpha \end{array}\right)}_{F}\cdot \left(\begin{array}{c}n\\y \end{array}\right)

F=[-beta 1;1 alpha]
%es igual a
lft(F_paralelo,F_serie,1,1)

% Para calcular la disipación "input" (resistencia mínima) u "output" (conductancia mínima) podemos comprobar que, a partir de de m=-\beta n+y, u=n+\alpha y tenemos:
% \beta u^2+\alpha m^2=(\beta n^2+\beta\alpha^2y^2+2\beta\alpha ny)+(\alpha\beta^2n^2+\alpha y^2-2\alpha\beta ny ) \\ =\beta n^2+\beta\alpha^2y^2+\alpha\beta^2n^2+\alpha y^2=(1+\alpha\beta)(\beta n^2+\alpha y^2)
% Dividiendo por 1+\alpha\beta:
% \frac{\beta u^2+\alpha m^2}{1+\alpha\beta}=\alpha y^2+\beta n^2
% Con lo cual:
% mn=(y-\beta n)n=y(u-\alpha y)-\beta n^2=yu-\alpha y^2-\beta n^2=yu-\frac{\beta u^2+\alpha m^2}{1+\alpha\beta}
% Con lo que, como \int mn\,dt\geq 0, tenemos:
% \int uy\,dt\geq \alpha \int y^2\,dt+\beta \int n^2\,dt=\frac{\beta}{1+\alpha\beta}\int u^2\,dt+\frac{\alpha}{1+\alpha\beta}\int m^2\,dt
% De modo que, ante una entrada fuente de intensidad, acotamos:
% \int uy\,dt\geq \alpha \int y^2\,dt = R_1 \int i^2\,dt
% *el lado de la derecha es cuando V_p=n=0, z es un cortocircuito
% Ante una entrada fuente de tensión, acotamos:
% \int uy\,dt\geq \frac{\beta}{1+\alpha\beta}\int u^2\,dt=\frac{1/R_2}{1+R_1/R_2}\int V^2\,dt=\frac{1}{R_1+R_2}\int V^2\,dt
% *el lado de la derecha es cuando i_p=m=0, z es un circuito abierto
%% Estabilidad robusta
% Pasemos a ver el lugar geométrico de la respuesta en frecuencia del sistema conectado al circuito para que el teorema de pasividad garantice estabilidad robusta:

figure
dibujatodo(F)
xline(0), xlim([-.1 .6])
hold on
plot(1/Rmin,0,'xr','LineWidth',3), plot(1/Rmax,0,'xk','LineWidth',3)
hold off

%% Criterio del círculo (realim +)
% Si la respuesta en frecuencia de G está fuera del círculo cuyo diámetro está entre 1/R_{\max} y 1/R_{\min},  y da tantas vueltas (antihorario) a ese círculo (bueno, criterio de  Nyquist sólo dice a +1/R_{min}, pero no puede tocar el círculo) como polos inestables de G en bucle abierto, entonces el bucle G-F-\Delta es estable.
% Recuperamos casos anteriores cuando \beta\to 0 ( R_2\to\infty, circuito serie, caso 1, cruz "negra" hacia cero) o \alpha\to0 ( R_1\to 0, circuito paralelo, caso 2, cruz "roja" hacia infinito).
% *En interconexión con realimentación negativa, el círculo sería en el lado izquierdo del plano complejo. También habría otras versiones del criterio del círculo si se hubiese considerado la entrada V y la salida i.
%
%% Conclusiones
% El criterio del círculo es un resultado famoso históricamente que generaliza el criterio de Nyquist (dar vueltas a un "punto" -1 o +1, según signo de realimentación) a dar vueltas a un círculo ante determinada configuración de pasividad expresable como LFT... Su analogía eléctrica es el circuito serie-paralelo; existe también una analogía mecánica.
% Como corolario, si z es una no-linealidad de primer y tercer cuadrante (estático), el resultado son no-linealidades con ganancia entre R_{min} y R_{max} (no-linealidad de sector).
%
% Nota: si R_{max}<\infty, el criterio del círculo también se puede demostrar a partir del teorema de pequeña ganancia si \beta>0, asumiendo como hipótesis de partida que \tilde \Delta=\frac{\Delta-R_{med}}{\rho} es contractivo, siendo R_{med}=\frac{1}{2}(R_{max}+R_{min}), \rho=\frac{1}{2}(R_{max}-R_{min}). Dicha derivación a partir de pequeña ganancia no entra en los objetivos de este material.
%
%
%% Funciones auxiliares

function dibujatodo(F)
disp('El valor de lft(F,G) es:')
H=lft(F,sym('G'))
disp('Por comprobar, inv(lft(1/H,F)) debe dar G:')
simplify(inv(lft(1/H,F)))
regrid=fliplr(logspace(-2.5,2,65));
imgrid=logspace(-3,2.5,65);
imgrid=[-fliplr(imgrid) 0 imgrid];
Nr=length(regrid);
Ir=length(imgrid);
disp('El semiplano izquierdo de H mapea a la zona verde')
for rek=1:Nr
    partereal=-regrid(rek);
    gr=[];
    for imk=1:Ir
        parteim=imgrid(imk);
        lG=1/lft(partereal+j*parteim,F);
        gr=[gr lG]; 
        if(rek==Nr)
            plot(gr,'Color',[0.1 .6 0.2],'LineWidth',4)
        else
            plot(gr,'g')
            if(rek==1)
                hold on
            end
        end
    end
end
hold off
axis equal
grid on
xlabel('Real'), ylabel('Imag')
end