%% Modelado de un circuito eléctrico (LFT) con DOS parámetros físicos inciertos.
% © 2023, Antonio Sala Piqueras. Universitat Politècnica de València. Todos los derechos reservados.
% Vídeo-presentación en https://youtu.be/jDPdPuJyeA0
%
%
% Parámetros nominales:

C1_0=1e-3;C2=2e-3; R1=8;R3=16; R4=32; R2_0=4;
syms dv1 dv2 V1 V2 incR2 incC1 Vm u i1m

%
%% standard model
% La resistencia "real" R2 es R2_0+incR2, el condensador "real" es C1_0+incC1 modelo es:

R2=R2_0+incR2;C1=C1_0+incC1;
 Model1=[C1*dv1==(u-V1)/R1-i1m; ...
   C2*dv2==(Vm-V2)/R3; ...
   i1m==(Vm-V2)/R3+Vm/R4; (V1-Vm)==R2*i1m ]
Unknowns=[dv1,i1m,Vm,dv2]; %remaining V1,V2,incR2,u
SolvedModel=solve(Model1,Unknowns)
StateEqs=simplify([SolvedModel.dv1; SolvedModel.dv2])
A=jacobian(StateEqs,[V1,V2])
B=jacobian(StateEqs,[u])

% Introduzcamos ecuación de salida... por ejemplo, se dispone de sensor/objetivos de control en V1 y Vm

OutputEqs=[V1;SolvedModel.Vm];
C=jacobian(OutputEqs,[V1 V2])
D=jacobian(OutputEqs,[u])

% El modelo "nominal" sería:

incR2=0;incC1=0;
Anum=eval(A)
Bnum=eval(B)
Cnum=eval(C)
Dnum=eval(D)
circuitnominal=ss(Anum,Bnum,Cnum,Dnum);

%% introducing changes in the model to incorporate LFT
% Redeclaremos incR2 como simbólica (borremos su valor de cero de antes), y introduzcamoslo en forma LFT

syms incR2 incC1 w_deltaR w_deltaC  Ana   %extra symbolic variables
Model_sinR2C1=[ dv1==Ana; ...
   C2*dv2==(Vm-V2)/R3; ...
   i1m==(Vm-V2)/R3+Vm/R4;];
EqResistor2_LFT= ( (V1-Vm)==R2_0*i1m+w_deltaR  );
EqC1_LFT= ( C1_0*Ana+w_deltaC==(u-V1)/R1-i1m );
Eq_Uncertainty = [w_deltaR==incR2*i1m; w_deltaC==incC1*Ana];

% La parte lineal "cierta" del modelo es esta:

Model2=[Model_sinR2C1; EqResistor2_LFT; EqC1_LFT]

% El modelo nominal será w_deltaR=w_deltaC=0.
% El modelo completo será (equiv. a Model 1):

Model3=[Model2; Eq_Uncertainty]
Unknowns3=[Unknowns w_deltaR w_deltaC Ana ];
SolvedModel3=solve(Model3, Unknowns3)

% Que resulta, en tanto a dv1, dv2, i1m, Vm, exactamente igual que el Model1. Comprobemos, por ejemplo A

LaMismaA=simplify(jacobian([SolvedModel3.dv1, SolvedModel3.dv2],[V1 V2]))
A-LaMismaA

%
%% Modelado LFT: NO usar Eq_Uncertainty, con lo que w_delta es una entrada "virtual"

SolvedModelLFT=solve(Model2,[Unknowns Ana])

% El resultado, claro, depende de w_deltaR y w_deltaC

StateEqsLFTModel=[SolvedModelLFT.dv1;SolvedModelLFT.dv2]

% Representación interna completa como upper LFT

A=jacobian(StateEqsLFTModel,[V1,V2])
GeneralisedInput=[w_deltaR, w_deltaC, u];
B_lft=jacobian(StateEqsLFTModel, GeneralisedInput)
GeneralisedOutputEq=[SolvedModelLFT.i1m;SolvedModelLFT.Ana;V1;SolvedModelLFT.Vm]; %añadimos una salida más al bloque incierto
C_lft=jacobian(GeneralisedOutputEq,[V1,V2])
D_lft=jacobian(GeneralisedOutputEq,GeneralisedInput)

% Comprobemos que la interconexión LFT da el resultado correcto (modelo original):

incR2=0; incC1=0;
sysuncert=diag([incR2;incC1]);
sys2=ss(eval(A),eval(B_lft),eval(C_lft),eval(D_lft));
circuito_otra_vez=lft(sysuncert,sys2);
norm(circuito_otra_vez-circuitnominal,'inf')

% Simulemos variaciones de resistencia y capacidades al azar:

for i=1:10
    incR2=8*(rand()-.5); incC1=8e-4*(rand()-.5); 
    sysuncert=diag([incR2;incC1]);
    circmodif=lft(sysuncert,sys2);
    bode(circmodif)
    hold on
    [incR2 incC1*1e4 norm(circmodif-circuitnominal,'inf')]
end
hold off, grid on

% Márgen de estabilidad ante incertidumbre no estructurada (R2 no lineal, con retardo, dinámica...)

mgst=1/norm(sys2(1:2,1:2),'inf')
    

% El resultado da un margen "pequeño"... problema: se ha mezclado incR2 con incC1, tienen unidades distintas... Necesaria "ponderación" (normalización) para que los márgenes tengan significado, mgst no se sabe las unidades que tiene.