Serie di Fourier: problema termine A0 dell'uscita

[904]

allora mi sono calcolato gli autovalori della matrice positiva e ho creato una variabile lameq che vale 1 se questi sono uguali e una notlameq che è il contrario poi ho creato la matrice M che nel posto 1,1 ha il primo autovalore nel posto 1,2 ha lameq+notlameq* il secondo autovalore nel posto 2,1 ha 1 e nel posto 2,2 ha notlameq poi mi sono calcolato X0 sostituendo a t lo zero dalle soluzioni negative di vC e iL poi in u0 gli ingressi positivi a tempo 0 in X0p ho scritto l'equazione di stato usando X0 come stato e U0 come ingresso , in vCp0 e iLp0 ho sostituito invece le soluzioni positive di vC e iL in 0 m in dvCp0 e diLp0 ho derivato la soluzione positiva di vC e iL e gli ho sostituito 0 poi ho ricavato Nc ed Ni usando i termini calcolati in precedenza e formando il vettore che al primo posto ha X0p- dvCp0 e al secondo posto X0- vCp0 e infine ho calcolato le due costanti facendo l'inverso di M e moltiplicandolo per Nc ed Ni

[RenzoDF]

...mi sono calcolato X0 sostituendo a t lo zero dalle soluzioni negative di vC e iL poi in u0 gli ingressi positivi a tempo 0 in X0p ho scritto l'equazione di stato usando X0 come stato e U0 come ingresso ...

Fin qui sembrerebbe Ok, ma quando poi scrivi ...

... in vCp0 e iLp0 ho sostituito invece le soluzioni positive di vC e iL in 0 ... in dvCp0 e diLp0 ho derivato la soluzione positiva di vC e iL e gli ho sostituito 0 ...

.... non concordo; anche se non ho capito bene cosa intendi con quel "sostituito 0" , VCp0 e iLp0 rappresentano le variabili di stato e di conseguenza non possono presentare discontinuità nel passaggio da t<0 a t>0, di conseguenza devono essere calcolate sulla soluzione stazionaria negativa, non positiva, mentre le loro derivate dovranno essere calcolate dalle soluzioni per iC e vL relative alla topologia della rete per t>0 ma con vC=vCp0 e iL=iLp0 .

Le righe di codice errate potrebbero quindi essere le seguenti

Codice: Seleziona tutto

vCp0=subs(sol_pos_vC,t,0);
iLp0=subs(sol_pos_iL,t,0);
dvCp0=subs(diff(sol_pos_vC),t,0)
diLp0=subs(diff(sol_pos_iL),t,0)

ma è solo una mia ipotesi, dopo una giornataccia demoralizzante di lotta con gli "squali" ... ed è tutta da verificare.

La soluzione permanente di vc(t) è corretta in quanto è direttamente verificabile con un semplice calcolo:

a) visto infatti che per t > 0 l'induttore porta ad un cortocircuito fra i suoi due morsetti, la tensione sul generatore pilotato deve per forza portarsi a zero in quanto la tensione di comando è una sua frazione ma di segno opposto, sarà quindi anch'esso equivalente ad un cortocircuito

b) ne segue che supponendo di trasformare il parallelo di j2(t) e R6 nell'equivalente serie secondo Thevenin, con generatore



e quindi



e resistenza equivalente pari a R6, con una impedenza totale di maglia



avremo che la vC particolare a regime sarà proprio





La parte da controllare è quindi proprio quella transitoria.

[904]

eh ma se fosse come dici tu allora pure il prof ha sbagliato nel suo esempio :

Codice: Seleziona tutto

clear all
close all
clc
syms R1 R2 R3 k
syms vs is vC iL
syms ia ib ic vsw isw
syms L C Vs0 Vs1 w t phi Is0
e1=-vs+R3*iL+vsw+R3*ia;
e2=-vsw-R1*isw+vC+R1*ia;
e3=-vC+R2*ic+k*ic;
%% generatori
vs_m=Vs0+Vs1*cos(w*t+phi);
is_m=Is0;
FVs1=Vs1/sqrt(2)*exp(i*phi);
%% open
sol_open=solve(e1,e2,e3,'ia,ic,vsw');
iA_op=simplify(sol_open.ia);
iC_op=simplify(sol_open.ic);
vSW_op=simplify(sol_open.vsw);
i_cap_op=simplify(iA_op-iC_op+is-isw);
v_ind_op=vSW_op;
dvC_dt_op=subs(i_cap_op/C,isw,0)
diL_dt_op=subs(v_ind_op/L,isw,0)
%% closed
sol_closed=solve(e1,e2,e3,'ia,ic,isw');
iA_cl=simplify(sol_closed.ia);
iC_cl=simplify(sol_closed.ic);
iSW_cl=simplify(sol_closed.isw);
i_cap_cl=simplify(iA_cl-iC_cl+is-iSW_cl);
v_ind_cl=vsw;
dvC_dt_cl=subs(i_cap_cl/C,vsw,0)
diL_dt_cl=subs(v_ind_cl/L,vsw,0)
%% soluzione a regime per t<0: interruttore chiuso
% la corrente dell'induttore è costante (diL_dt_cl=0)
% la corrente is è costante
% si determina la tensione stazionaria del conensatore
% mediante l'ed. diff. dvC_dt_cl
vCp_m=simplify(subs(solve(dvC_dt_cl,'vC'),is,is_m));
isw_m_staz=solve(subs(iSW_cl,[vC,vsw,iL,vs],[vCp_m,0,isw,Vs0])-isw,'isw');
iLp_m=isw_m_staz;
pretty(iLp_m)
pretty(vCp_m)
%% condizioni iniziali
iL_0=iLp_m;
vC_0=vCp_m;
%% soluzione particolare per t>0
% componente stazionaria
sol_p_staz=solve(dvC_dt_op,diL_dt_op,'vC,iL')
vCp_p_staz=simplify(subs(sol_p_staz.vC,[is,vs],[is_m,Vs0]));
iLp_p_staz=simplify(subs(sol_p_staz.iL,[is,vs],[is_m,Vs0]));
pretty(iLp_p_staz)
pretty(vCp_p_staz)
% componente sinusoidale
sol_p_sin=solve(dvC_dt_op-i*w*vC,diL_dt_op-i*w*iL,'vC,iL')
FvCp_p_sin=simplify(subs(sol_p_sin.vC,[is,vs],[0,FVs1]));
FiLp_p_sin=simplify(subs(sol_p_sin.iL,[is,vs],[0,FVs1]));
pretty(FiLp_p_sin)
pretty(FvCp_p_sin)
%% risultati numerici
R1=1; R2=2; R3=2; k=2; C=50e-6; L=20e-3;
Is0=10; Vs0=5; Vs1=10;
w=1/sqrt(L*C); phi=pi/4; T=2*pi/w;
FvCp_p=subs(FvCp_p_sin)
FiLp_p=subs(FiLp_p_sin)
vCp_p_sin=abs(FvCp_p)*sqrt(2)*cos(w*t+angle(FvCp_p))
iLp_p_sin=abs(FiLp_p)*sqrt(2)*cos(w*t+angle(FiLp_p))
%% soluzioni a regime per t>0 complete
vCp_p=vCp_p_sin+subs(vCp_p_staz)
iLp_p=iLp_p_sin+subs(iLp_p_staz)
vLp_p=L*diff(iLp_p,'t')
% h1=figure(1)
% subplot(211), ezplot(vCp_p,[0 3*T])
% hold on, ezplot(vLp_p,[0 3*T]), axis auto, title('vCp,vLp t>0')
% subplot(212), ezplot(iLp_p,[0 3*T]), title('iLp t>0')
% dock
%% matrice dinamica A t>0
dvC=subs(dvC_dt_op);
diL=subs(diL_dt_op);
A(1,1)=eval(subs(dvC,[is,vs,iL],[0,0,0])/vC);
A(1,2)=eval(subs(dvC,[is,vs,vC],[0,0,0])/iL);
A(2,1)=eval(subs(diL,[is,vs,iL],[0,0,0])/vC);
A(2,2)=eval(subs(diL,[is,vs,vC],[0,0,0])/iL);
lam=eig(A)
%% matrice B t>0
B(1,1)=eval(subs(dvC,[iL,vC,is],[0,0,0])/vs);
B(1,2)=eval(subs(dvC,[iL,vC,vs],[0,0,0])/is);
B(2,1)=eval(subs(diL,[iL,vC,is],[0,0,0])/vs);
B(2,2)=eval(subs(diL,[iL,vC,vs],[0,0,0])/is);
%% condizioni iniziali
X0=subs([vC_0;iL_0])
U0=subs([is_m;subs(vs_m,'t',0)])
Xp0=A*X0+B*U0
vCpo_p=subs(vCp_p,'t',0);
iLpo_p=subs(iLp_p,'t',0);
Xpo_p=[vCpo_p;iLpo_p];
dvCp_pdt=diff(vCp_p,'t');
dvCpo_pdt=subs(dvCp_pdt,'t',0)
diLp_pdt=diff(iLp_p,'t');
diLpo_pdt=subs(diLp_pdt,'t',0)
%% calcolo coefficienti
Xo_p=X0;
dXdto_p=Xp0;
lameq=(lam(1)==lam(2))
nlameq=not(lameq)
M=[lam(1) lameq+lam(2)*nlameq;
   1      nlameq]
Nc=[dXdto_p(1)-dvCpo_pdt;Xo_p(1)-vCpo_p]
Ni=[dXdto_p(2)-diLpo_pdt;Xo_p(2)-iLpo_p]
kC=inv(M)*Nc
kI=inv(M)*Ni
%% costanti di tempo
tau1=-1/real(lam(1))
tau2=-1/real(lam(2))
%% soluzione completa
vC_t=kC(1)*exp(lam(1)*t)+kC(2)*exp(lam(2)*t)+vCp_p
iL_t=kI(1)*exp(lam(1)*t)+kI(2)*exp(lam(2)*t)+iLp_p
%% grafici delle soluzioni complete
% h2=figure(2);
% T1=T
% Tt=5*max(T1,max(tau1,tau2));
% subplot(211), ezplot(subs(vCp_m)*t/t,[-T1,0])
% hold on, ezplot(vC_t,[0,Tt]), axis auto, grid on
% subplot(212), ezplot(subs(iLp_m)*t/t,[-T1,0])
% hold on, ezplot(iL_t,[0,Tt]), axis auto, grid on
% set(gcf,'WindowStyle','Docked')
%% analisi di Fourier
% segnale di ingresso: onda quadra
syms AH AL
App=2; % valore picco-picco
Aav=1; % valore medio
D=0.25; % duty-cycle
fs=120; % frequenza fondamentale
ws=fs*2*pi; % pulsazione fondamentale
Ts=1/fs; % periodo
tH=D*Ts; tL=(1-D)*Ts;
% calcolo di AH ed AL
eq1=AH-AL-2;
%eq1='AH-AL=2'
eq2=AH*D+AL*(1-D)-1;
AA=solve(eq1,eq2,'AH,AL');
AH=(AA.AH)
AL=(AA.AL)
% calcolo dei coefficienti della serie di Fourier
Na=10; % numero di armoniche
ao_vs=eval(1/Ts*(int(AH,t,0,tH)+int(AL,t,tH,Ts)))
ar_vs=zeros(1,Na);
br_vs=zeros(1,Na);
vs_t=ao_vs;
h=WAITBAR(0,'..running..');
IA=eye(2,2);
ao_X=-inv(A)*B*[ao_vs;0];
ao_vc=ao_X(1);
vc_t=ao_vc;
for r=1:Na
    waitbar(r/Na);
    ar_vs(r)=2/Ts*(int(AH*cos(r*ws*t),t,0,tH)+int(AL*cos(r*ws*t),t,tH,Ts));
    br_vs(r)=2/Ts*(int(AH*sin(r*ws*t),t,0,tH)+int(AL*sin(r*ws*t),t,tH,Ts));
    temp=(ar_vs(r)+i*br_vs(r));
    cr_vs(r)=abs(temp);
    phir_vs(r)=-angle(temp);
    FVsr=cr_vs(r)*exp(i*phir_vs(r));
    FXr=inv(i*r*ws*IA-A)*B*[FVsr;0];
    FVcr=FXr(1);
    vs_t=vs_t+cr_vs(r)*cos(r*ws*t+phir_vs(r));
    vc_t=vc_t+abs(FVcr)*cos(r*ws*t+angle(FVcr));
end
figure, dock
ezplot(vs_t,[0,2*Ts]);
hold on
ezplot(vc_t,[0,2*Ts]);
axis auto
% ic_t=diff(vc_t);
% ezplot(ic_t,[0,2*Ts]);
% axis auto







[RenzoDF]

eh ma se fosse come dici tu allora pure il prof ha sbagliato nel suo esempio ...

A quale problema (rete) si riferisce questo esempio? ... direi comunque che non possiamo perdere tempo andando a controllare cosa ha scritto il tuo professore, siamo già troppo "incartati" con il tuo di codice

Tanto per curiosità, hai provato a simulare con le modifiche che ti ho suggerito? ... io come sai non posso farlo.

Se mi puoi postare i risultati numerici forniti da Matlab relativi alle seguenti righe del tuo codice :

Codice: Seleziona tutto

%transitorio
lam=eig(Apos);
lameq=(lam(1)==lam(2));
notlameq=not(lameq);
M=[lam(1) lameq+ notlameq*lam(2);1 notlameq];
X0=[subs(sol_neg_vC,t,0);subs(sol_neg_iL,t,0)];
U0=[v1pos;subs(j2,t,0)];
X0p=Apos*X0+Bpos*U0;
vCp0=subs(sol_pos_vC,t,0);
iLp0=subs(sol_pos_iL,t,0);
dvCp0=subs(diff(sol_pos_vC),t,0)
diLp0=subs(diff(sol_pos_iL),t,0)
Nc=[X0p(1)-dvCp0;X0(1)-vCp0];
Ni=[X0p(2)-diLp0;X0(2)-iLp0];
Kc=inv(M)*Nc;
Ki=inv(M)*Ni;
vC_t=sol_pos_vC+ Kc(1) * exp(lam(1)*t)+ Kc(1)*exp(lam(2)*t);
iL_t=sol_pos_iL+ Ki(1) * exp(lam(1)*t)+ Ki(1)*exp(lam(2)*t);
%%trovo le tensioni e traccio il grafico
vLneg=subs(sol_neg.vL,[vC iL vSw],[sol_neg_vC sol_neg_iL 0]);
vLpos=subs(sol_pos.vL,[v1_s iL iSw],[v1pos iL_t 0]);


ovvero M,X0, U0, X0p, vCp0, iLp0, dvCp0, diLp0, ...., vC_t e iL_t, mi faresti un favore.


BTW ... sarei poi anche curioso di capire il significato di queste due righe del tuo codice

Codice: Seleziona tutto

num=[75 0];
den=2*conv([2 1625],[1 50]);

ovvero, da dove arrivano quei valori numerici.

[904]

appena posso ti posto l output comunque quei due den e num arrivano dall output sono il numeratore e denominatore della funzione di rete la traccia relativa a quell esempio del prof non la so perché il prof ha messo gli appunti senza la traccia relativa a quell esercizio ma da quello che spiegò la parte del transitorio è standard è sempre quello il transitorio non si modifica quel procedimento dice lui

[904]

Ho capito la parte sul transitorio andando a ricevimento dal professore
Allora vi posto una diapositiva che potrà aiutarvi a capire anche a voi perché avevi ragione tu ci vuole la soluzione negativa

wwww.png (647.63 KiB) Osservato 2641 volte

BTW: il sito da chrome non mi va solo da firefox

[RenzoDF]

Ho capito la parte sul transitorio andando a ricevimento dal professore
Allora vi posto una diapositiva che potrà aiutarvi a capire anche a voi perché avevi ragione tu ci vuole la soluzione negativa

Ovvero che ti ha detto il professore?

[904]

mi ha spiegato che per calcolare il transitorio basta che noi nella soluzione del equazione differenziale che sarà fatta da una soluzione particolare che è quella a regime e da un'altra che dipende dal transitorio è sarà nella forma come quella di tutte le equazioni omogenee da qui per calcolare le costanti si porta al primo membro la soluzione particolare e si impongono le condizioni iniziali e cioè si calcola per esempio vC in 0- e la sua derivata in 0+ mentre la soluzione particolare in 0+ non ho capito però il perché di questo vC in 0- la sua derivata in 0+ e la soluzione particolare in 0+ lui ha detto che vC lo calcola in 0- perché dipende dallo stato precedente del circuito la matrice M è una generalizzazione che ha fatto lui e non dobbiamo usarla sennò rischiamo di non capire il transitorio

[RenzoDF]

Io non vorrei essere ripetitivo, ma se quello che ti ha detto corrisponde a quello che hai scritto, non ti ha detto nulla di interessante.

Ripeto: le condizioni iniziali dobbiamo calcolarcele entrambe a t=0+

a) per la condizione iniziale di vC, si deve usare la vC(0+); e solo in quanto vC NON può presentare discontinuità si usa vC(0+)=vC(0-) e fin qui non ci piove ... e lo ha detto pure lui, ma la parte più importante è la seguente ...

b) per la condizione iniziale sulla derivata v'C(0+) non si può assumere allo stesso modo che sia uguale a v'C(0-) in quanto iC non presenta il vincolo di continuità, ne segue che la v'C(0+) dobbiamo ricavarcela in un altro modo, ovvero attraverso il calcolo della iC(0+) in una rete che rappresenti la situazione istantanea al tempo t=0+.

Supporremo cioè di "congelare la rete a t=0+", sostituendo a C un GIT di valore pari a vC(0+) e a L un GIC di valore pari a iL(0+) e ricavare iC(0+) e vL(0+), che rappresentano a meno di una costante 1/C ed 1/L rispettivamente, le derivate della tensione e della corrente all'istante del congelamento t=0+.

[904]

si ma il problema è perché dobbiamo calcolarcele a t=0+?