%% Pasividad + LFT: ampliación modelado
% © 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/pasexc2.html
% Este código funcionó con Matlab R2020a
% Objetivo: Ampliar el modelado de interconexiones LFT con otros escenarios para aplicar pasividad, y discutir analogía con ejemplos mecánicos.
% Tabla de Contenidos
% Revisión teórica
% Impedancias serie o paralelo
%     Caso 1:  exceso de pasividad input
%     Caso 2, exceso de pasividad output
%     Unificación caso 1 y caso 2
% 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, 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

%
%
%% 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 1 -2.5 2.5])

% G debe ser estable.
%% Unificación caso 1 y caso 2
% Si se eligiera V como entrada e i como salida (o F como entrada y v como salida en el caso mecánico), entonces los sistemas eléctricos o mecánicos del caso 1 servirían como ejemplo del caso 2 y viceversa.
%
% El teorema de pasividad (región del plano complejo sigue significando lo mismo...) en el caso 1 (izquierda), G debe calcular una intensidad a inyectar en función de medidas de tensión; en el caso 2 (derecha), G debe calcular una tensión a aplicar en función de medidas de intensidad. Dependiendo del escenario técnico será de aplicación uno u otro (o ninguno de ellos estrictamente, siendo parte de una planta generalizada mayor).
%
%%
%%
%
%% 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,0,25));
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