27/11/2013

%% Prova d'esame del 27/11/2013
%% Alise Dario 0612701333

%% Operazioni di pulizia

clc
clear all
close all

%% Dati

syms k t iL

R1 = 1; R2 = 1; R3 = 1;
micro = 1e-6;
C = 10*micro;
L = 10*micro;

%% Generatori

JS2 = k*iL;
amp = 10;
wg = 100;
phi = pi/4;
VS = amp*cos(wg*t+phi);

%% Analisi per t<0

% Meotodo dei potenziali nodali

syms Ua vC vL iC iL Vs

sol = solve((Ua+vC)/R3-iL+(Ua-Vs)/R1, 'Ua');
Ua = sol;
vL_m = -Ua;
iC_m = solve((vC+Ua)/R3+(Vs+vC)/R2+iC, 'iC');
dvCdt_m = iC_m/C;
diLdt_m = vL_m/L;

% Modello di stato

% Matrice A

A11 = subs(dvCdt_m, [vC, iL, Vs], [1, 0, 0]);
A12 = subs(dvCdt_m, [vC, iL, Vs], [0, 1, 0]);
A21 = subs(diLdt_m, [vC, iL, Vs], [1, 0, 0]);
A22 = subs(diLdt_m, [vC, iL, Vs], [0, 1, 0]);

A_m = [A11, A12; A21, A22];

% Vettore B

B1 = subs(dvCdt_m, [vC, iL, Vs], [0, 0, 1]);
B2 = subs(diLdt_m, [vC, iL, Vs], [0, 0, 1]);

B_m = [B1; B2];

% Autovalori per la stabilità

lam_m = eig(A_m)

% Calcolo delle soluzioni di regime (unicamente soluzione sinusoidale
% ricavata con il metodo dei fasori)

Id = eye(2);
FVS = amp*exp(j*phi);
FX = inv(j*wg*Id-A_m)*B_m*FVS;

vC_m = abs(FX(1))*cos(wg*t+angle(FX(1)));
iL_m = abs(FX(2))*cos(wg*t+angle(FX(2)));
X_m = [vC_m; iL_m];
X_m_0 = subs(X_m, 't', 0);

% Calcolo della grandezza di interesse

IR3_m = subs((Ua+vC)/R3, [vC, iL, Vs], [vC_m, iL_m, VS]);
PR3_m = (IR3_m)^2*R3;

%% Analisi per t>0

clear Ua vC vL iC iL Vs
syms Ua Ub vC vL iC iL Vs

eq1 = (Ua-Vs)/R1-iL+(Ua+vC)/R3;
eq2 = -k*iL+(Ub+vC)/R2;
sol = solve(eq1, eq2, 'Ua, Ub');
vL_p = -sol.Ua;
eq1 = (sol.Ua+vC)/R3+k*iL+iC;
iC_p = solve(eq1, 'iC');
dvCdt_p = iC_p/C;
diLdt_p = vL_p/L;

% Modello di stato

% Matrice A

A11 = subs(dvCdt_p, [vC, iL, Vs], [1, 0, 0]);
A12 = subs(dvCdt_p, [vC, iL, Vs], [0, 1, 0]);
A21 = subs(diLdt_p, [vC, iL, Vs], [1, 0, 0]);
A22 = subs(diLdt_p, [vC, iL, Vs], [0, 1, 0]);

A_p = [A11, A12; A21, A22];

% Vettore B

B1 = subs(dvCdt_p, [vC, iL, Vs], [0, 0, 1]);
B2 = subs(diLdt_p, [vC, iL, Vs], [0, 0, 1]);

B_p = [B1; B2];

% Autovalori per la stabilità

lam_p = eig(A_p)

% Per la stabilità k deve essere k>-1/2
% Pongo k = 1;
k = 1;
A_p = subs(A_p);
B_p = subs(B_p);
lam_p = subs(lam_p);

% Calcolo delle soluzioni particolari

FVS = amp*exp(j*phi);
FX = inv(j*wg*Id-A_p)*B_p*FVS;
vCP_p = abs(FX(1))*cos(wg*t+angle(FX(1)));
iLP_p = abs(FX(2))*cos(wg*t+angle(FX(2)));

% Soluzione generale

syms KC1 KC2 KL1 KL2

vCG_p = KC1*exp(lam_p(1)*t)+KC2*exp(lam_p(2)*t);
iLG_p = KL1*exp(lam_p(1)*t)+KL2*exp(lam_p(2)*t);

% Soluzione completa

vC_p = vCG_p + vCP_p;
iL_p = iLG_p + iLP_p;

% Calcolo delle costanti

U = VS;
U_0 = subs(U, 't', 0);
dX0_dt = A_p*X_m_0+B_p*U_0;

KC1 = solve(subs(vC_p, 't', 0)-X_m_0(1), 'KC1');
vC_p = subs(vC_p);
KC2 = solve(subs(diff(vC_p), 't', 0)-dX0_dt(1), 'KC2');
vC_p = subs(vC_p);

KL1 = solve(subs(iL_p, 't', 0)-X_m_0(2), 'KL1');
iL_p = subs(iL_p);
KL2 = solve(subs(diff(iL_p), 't', 0)-dX0_dt(2), 'KL2');
iL_p = subs(iL_p);

% Calcolo della grandezza di interesse

IR3_p = subs((sol.Ua+vC)/R3, [vC, iL, Vs], [vC_p, iL_p, VS]);
PR3_p = (IR3_p)^2*R3;

% Grafici

max = max(eval(lam_p));

figure(1)
subplot(211), ezplot(PR3_m, [-5*pi/wg, 0]), grid on
title('Potenza istantanea resistore R3 con t<0')
subplot(212), ezplot(PR3_p, [0, 5*pi/wg]), grid on
title('Potenza istantanea resistore R3 con t>0')

%% Frequenza di risonanza

s = tf('s');

% Caso t<0

Hs_m = inv(s*Id-eval(A_m))*eval(B_m);
HCs_m = zpk(minreal(Hs_m(1)));
figure(2), bode(HCs_m), grid on

syms w
G = -150000; Z = 3.333e004; P = 1e005;
MHw_m = abs(G)*sqrt((w^2+Z^2)/(w^2+P^2)^2);
MHw_m2 = MHw_m^2;
wr = eval(solve(diff(MHw_m2), 'w'));
wr_m = wr(2)

% Caso t>0

Hs_p = inv(s*Id-eval(A_p))*eval(B_p);
HCs_p = zpk(minreal(Hs_p(1)));
figure(3), bode(HCs_p), grid on

G = -50000; Z = -1e005; W02 = 1e010; R = 1e005;

MHw_p = abs(G)*sqrt((w^2+Z^2)/((W02-w^2)^2+w^2*R^2));
MHw_p2 = MHw_p^2;
wr = eval(solve(diff(MHw_p2), 'w'));
wr_p = wr(2)

%% Potenza massima pilotato

IP = k*iLP_p;
VP = R2*IP-vCP_p;
P = IP*VP;
figure(4), ezplot(P,[0,5*pi/wg]), grid on, title('Potenza pilotato t>0')

% Considero esclusivamente la soluzione particolare, da cui si evidenzia un
% massimo

tmax = eval(solve(diff(P,'t')));
t0 = tmax(2)

% La potenza massima è all'istante t0, per cui sostituisco

digits(3)
Pmax = eval(subs(P, 't', t0))

%% Percentuale componente media

VM = 0.25;
X0 = inv(A_p)*(-B_p)*VM;
VC0 = X0(1);

FF = 500;
T = 1/FF;
wq = 2*pi/T;
AH = 1;
AL = 0;
dc = 0.25;
TH = T*dc;
perc = 0.05;

syms n

An = (2/T)*int(AH*cos(n*wq*t), t, 0, TH);
Bn = (2/T)*int(AH*sin(n*wq*t), t, 0, TH);
FCn = An-j*Bn;

HXn_p = inv(j*wq*n*Id-A_p)*(B_p);
Hn_p = HXn_p(1);
H0_p = abs(subs(Hn_p,'n',0));

Hlim = perc*H0_p;
wlim = solve(MHw_p2-Hlim^2, 'w');
nmax = round(wlim(1)/wq);

for n=1:nmax
    FVC = subs(Hn_p)*subs(FCn);
    if(abs(FVC)>perc*VC0)
    nmax = n;
    end
end

% Il mio ragionamento è questo: poichè ho determinato l'nmax che
% rappresenta il numero massimo di armoniche del generatore in ingresso che
% non sono inferiori al 5% della sua componente media sono sicuro che se
% controllando le ampiezze dell'uscita (la tensione del condensatore) salvo
% l'indice dell'ultima ampiezza avente modulo maggiore della percentuale di
% componente media di interesse mi risparmio tutte quelle osservazioni
% grafiche che possono risultare insidiose

VCQ = VC0;
for n=1:nmax
    FVCn = subs(Hn_p)*subs(FCn);
    VCQ = VCQ+abs(FVCn)*cos(n*wq*t+angle(FVCn));
end

% Il grafico della tensione del condensatore costituito dalle
% prime 46 armoniche è riportato in Figura E9
figure(5), ezplot(VCQ,[0,2*T]), grid on, title('vC(t)')