%% Modelado circuito eléctrico
% (c) 2022 A. Sala, Universitat Politècnica de Valencia. Todos los derechos reservados.
% Vídeo-presentación disponible en: http://personales.upv.es/asala/YT/V/cir1.html
%
%
% declaremos los símbolos a usar:

syms U1 U2 VR1 iR1 dVC1 VC1 iC1 dVC2 VC2 iC2 diL VL VR2 iL VR3 iU2 real

% y los parámetros constantes:

R1=100; R2=100; R3=100; C1=1e-3;C2=2e-3; L=2e-2;

% Introduzcamos las ecuaciones de la física (Nota: La sintaxis es diferente en otras versiones del Symbolic Toolbox; en este caso se ha usado Matlab R2016a):

dinamica=[dVC1==iC1/C1; dVC2==iC2/C2; diL==VL/L];
resistencias=[VR1==iR1*R1;  VR2==iL*R2; VR3==iC1*R3];
balances=[iR1==iU2+iC1; iU2==iL+iC2; U1==VC1+VR3+VR1; VR3+VC1==VC2-U2; VL+VR2==VC2];
Modelo =[dinamica; resistencias ; balances];
Modelo

% Comprobemos cuántas ecuaciones tiene el modelo:

length(Modelo)

% y cuántas variables simbólicas aparecen:

symvar(Modelo)
length(symvar(Modelo))

% Se pueden "despejar", por tanto, 11 símbolos en función de otros cinco.
% Para obtener la representación interna (EDO) normalizada:

 Incognitas=[iR1, VR1, VR2, VR3, diL, VL, iC1, dVC1, iC2, dVC2, iU2];
 length(Incognitas)

% Podemos resolver el sistema de 11 ecs con 11 incognitas, en funcion de las cinco letras restantes  que son (U1, U2, VC1, VC2, iL)

 sol=solve(Modelo, Incognitas);
 

% Como no ha dado ningun error, el modelo, era, en efecto, un sistema de ecuaciones algebraico-diferencial de índice 1 dado que hemos podido despejar todo en función de las variables de estado y las entradas:

%Entradas
VectorEntradas=[U1;U2];
 % vbles estado, aquellas bajo signo derivada:
VectorEstado=[VC1;VC2;iL];

% Las ecuaciones de estado normalizadas son:

ecs_estado=[sol.dVC1;sol.dVC2;sol.diL]

% Las ecuaciones de salida son cualesquiera de las incógnitas de interés... por ejemplo:

%solo nos interesa VC2, la intens por R1, y la intens. por fuente U2
ecs_salida=[VC2; sol.iR1;sol.iU2]

% Nota: VC2 no necesita "despejarse" en funcion del estado y entradas porque directamente ES una vble de estado. (poner sol.VC2 sería un error, porque NO es incógnita).
% En forma matricial:

A=jacobian(ecs_estado,VectorEstado)
B=jacobian(ecs_estado,VectorEntradas)
C=jacobian(ecs_salida,VectorEstado)
D=jacobian(ecs_salida,VectorEntradas)

% Con lo que podríamos expresar el modelo normalizado:
% \dot x=Ax+Bu \\
% y=Cx+Du
% Siendo x un vector de longitud 3 (VectorEstado), u de longitud 2 (VectorEntradas), e y de longitud 3 (y=(VC2, iR1,iU2)').
% Vamos a simular, por ejemplo, el transitorio con alimentación CC (U1, U2 constantes):

modelocircuito=ss(eval(A),eval(B),eval(C),eval(D));
% eval es para pasar de simbólico (cadena de caracteres/texto) a numérico
modelocircuito.OutputName={'VC2','iR1','iU2'};
modelocircuito.InputName={'U1','U2'};
step(modelocircuito)