%% Modelování jednoduchého RLC obvodu: stavový model, přenosová funkce, algebro-diferenciální rovnice

%%
% Uvažujme jednoduchý RLC obvod na obrázku
%
% <<dae_circuit.png>>
%

%% Parametry
l = 0.0015;     % [H]
c = 1e-6;       % [F]
r1 = 100;       % [Ohm]
r2 = 20;        % [Ohm]
u0 = 10;        % [V]

%% Sestavení matic stavového modelu
A = [-(r1+r2)/(r1*r2*c), 0; 0, 0];
B = [1/(r1*c); 1/l];
C = [1/r2, 0];
D = 0;

%% Vytvoření stavového modelu
Gss = ss(A,B,C,D);

set(Gss,'InputName','u_o')
set(Gss,'OutputName','i_2')
set(Gss,'StateName',{'u_C','i_L'})


%% Analýza dynamiky výpočtem vlastních čísel
damp(Gss)

%% Simulace odezvy na změnu vstupního napětí
step(Gss)

%% Simulace odezvy na změnu vstupního napětí voláním obecného solveru pro nelineární ODE
t1 = 1.5e-4;
tspan = [0 t1];
u = 1;
x0 = [0;0];
options = odeset('RelTol',1e-4,'AbsTol',1e-6);

[t,x] = ode45(@(t,x) A*x+B*u,tspan,x0,options);
y = x(:,1)/r2;

plot(t,y)
xlabel('t')
legend({'i_2'})

%% Přenosová funkce pro obvod
Gtf = tf(Gss)

step(Gtf)

%% 
% Nedá se nevšimonut si, že řád přenosové funkce je 1, zatímco stavový
% model je řádu dva (jsou v něm dva dynamické prvky - cívka a kondenzátor).
% Jak je to možné?

tzero(Gss)

%% Vytvoření matic pro lineární deskriptorový (také implicitní DAE) model
E = zeros(9,9);
E(1,1) = l;
E(2,2) = c;

A = zeros(9,9);
A(1,9) = 1;
A(2,8) = 1;
A(3,4) = -r1; A(3,5) = 1;
A(4,6) = -r2; A(4,7) = 1;
A(5,2) = -1; A(5,5) = -1; 
A(6,5) = 1; A(6,7) = 1; A(6,9) = -1;
A(7,2) = 1; A(7,7) = -1;
A(8,1) = -1; A(8,3) = 1; A(8,4) = -1;
A(9,4) = 1; A(9,6) = -1; A(9,8) = -1;

B = zeros(9,1);
B(5,1) = 1;

C = zeros(1,9);
C(1,6) = 1;

D = 0;

%% Sestavení deskriptorového modelu v Matlabu
Gdae = dss(A,B,C,D,E);
set(Gdae,'InputName','u_o')
set(Gdae,'OutputName','i_2')
set(Gdae,'StateName',{'i_L','u_c','i_0','i_1','u_1','i_2','u_2','i_c','u_L'})

%% Analýza deskriptorového modelu
damp(Gdae)

%% Simulace odezvy na jednotkový skok ve vstupním napětí
step(Gdae)

%% Simulace odezvy na změnu vstupního napětí voláním obecného solveru pro nelineární ODE
t1 = 1.5e-4;
tspan = [0 t1];
u = 1;
z0 = zeros(9,1);
M = zeros(9,9); M(1,1) = l; M(2,2) = c;
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);

[t,z] = ode15s(@(t,z) A*z+B*u,tspan,z0,options);
y = z(:,6);

plot(t,y)
xlabel('t')
legend({'i_2'})

%% Převod deskriptorového modelu na stavový model

E11 = E(1:2,1:2);
E12 = E(1:2,3:end);
E21 = E(3:end,1:2);
E22 = E(3:end,3:end);

A11 = A(1:2,1:2);
A12 = A(1:2,3:end);
A21 = A(3:end,1:2);
A22 = A(3:end,3:end)

B1 = B(1:2,:);
B2 = B(3:end,:);

%%
% A s tímto novým (permutovaným) stavovým popisem můžu eliminovat (pokud to
% půjde) ty algebraické proměnné a ponechat v modelu jen ty diferenciální

As = E11\A11 - E11\A12*(A22\A21)
Bs = E11\B1 - E11\A12*(A22\B2)

%%  Analýza stavového modelu získaného naším vlastním převodem z DAE
eig(As)

%% Modifikace obvodu - prohození L a C
% Uvažujme teď jiný obvod - u toho předchozího prohodíme prvky L a C
%
% <<dae_circuit_singular.png>>

A = zeros(9,9);
A(1,9) = 1;
A(2,8) = 1;
A(3,4) = -r1; A(3,5) = 1;
A(4,6) = -r2; A(4,7) = 1;
A(5,7) = -1; A(5,9) = -1; 
A(6,2) = 1; A(6,5) = -1; A(6,7) = -1;
A(7,7) = 1; A(7,9) = -1;
A(8,3) = 1; A(8,4) = -1; A(8,9) = -1;
A(9,1) = -1; A(9,4) = 1; A(9,6) = -1;

%% Sestavení deskriptorového modelu v Matlabu
Gdae = dss(A,B,C,D,E);
set(Gdae,'InputName','u_o')
set(Gdae,'OutputName','i_2')
set(Gdae,'StateName',{'i_L','u_c','i_0','i_1','u_1','i_2','u_2','i_c','u_L'})

%% Analýza deskriptorového modelu
damp(Gdae)

%% Simulace odezvy na jednotkový skok ve vstupním napětí
step(Gdae)

%% Převod deskriptorového modelu na stavový model
A11 = A(1:2,1:2);
A12 = A(1:2,3:end);
A21 = A(3:end,1:2);
A22 = A(3:end,3:end)

%%
% V této fázi by bylo na místě začít vyjadřovat těch 7 semi-stavových proměnných 
% coby funkce těch pouhých dvou stavových proměnných i_L a u_C. Jenže tento postup selhává, 
% proto odpovídající matice je singulární. 
rank(A22)

%%
% Zjevně jde o problém indexu vyššího než 1. To nakonec půjde vidět, když
% zkusíme ten problém řešit numericky pomocí obecného solveru pro (tuhé)
% obyčejné diferenciální rovnice.

%% Simulace odezvy na změnu vstupního napětí voláním obecného solveru pro nelineární ODE
t1 = 1.5e-4;
tspan = [0 t1];
u = 1;
z0 = zeros(9,1);
M = zeros(9,9); M(1,1) = l; M(2,2) = c;
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);

[t,z] = ode15s(@(t,z) A*z+B*u,tspan,z0,options);
y = z(:,6);

plot(t,y)
xlabel('t')
legend({'i_2'})