% DEFINIÇÃO DE PARAMETROS DE SIMULAÇAO clear close all fr=50; % FRECUENCIA DESEJADA DA TENSAO DE SAIDA EM Hz w=round(2*pi*fr); % FRECUENCIA ANGULAR DESEJADA DA TENSAO DE SAIDA EM rad/seg mf=160; % NUMERO DE PULSOS PWM POR PERIODO fs=mf*fr; % FREQUENCIA DE FREQUENCIA DE COMUTAÇÃO Ts=1/fs; % PERIODO DE COMUTAÇÃO) fa=1*fs; % FREQUENCIA DE AMOSTRAGEM Ta=1/fa; % PERIODO DE AMOSTRAGEM e atualizaçao da malha de precompensaçao da planta. frep=fs; % Frequencia de amostragem e ftualiazaçao da malha de Controle com Modelo Interno Tr=1/frep; % PERIODO Atualiazaçao da Açao de Controle do Modelo Interno ms=round(fa/fr); % NUMERO DE AMOSTRAS POR PERIODO DA TENSAO DE SAIDA mr=round(frep/fr); % NUMERO DE AMOSTRAS DO cONROLADOR DE MODELO INTERNO EM PERIODO DA TENSAO DE SAIDA nc=200; % NUMERO DE PERIODOS DE SIMULAÇAO 200 na=nc*ms; % NUMERO DE ToTAL DE AMOSTRAS POR SIMULAÇAO p=1024; % NUMERO DE PONTOS POR AMOSTRA (PASSO DE SIMULAÇAO) deve ser um numenro par np=na*p; % NUMERO TOTAL DE PONTOS DE SIMULAÇAO nr=mr*nc; % NUMERO TOTAL DE amostras do controlador do modelo interno em uma simulaçao Vbase=311; % VALOR BASE DE TENSAO PARA NORMALIZAÇAO Ibase=50; % VALOR BASE DE CORRENTE PARA NORMALIZAÇAO Vcc=400; N=mr; % Numero de amostras do controlador do modelo interno em um periodo da tensao de saida % e tambem o ordem do controlador do modelo interno % sem incluir o filtro Q e um possil compensador. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DA PLANTA. Filtro e Carga Lf=1e-3; % Indutancia de entrada em uH Cf=60e-6; % Capacitancia do filtro de saida em uF wo=1/sqrt(Lf*Cf); % Frequencia angular de resonancia do filtro LC em rad/seg fo=wo/2/pi; % Frequencia de resonancia do filtro LC em Hz Rc=10; % Resistencia de carga %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% s=tf('s'); L=Lf; C=Cf; R=Rc; Gpc=(1/(L*C))/(s^2 + s/(R*C) + 1/(L*C)); wn=1/sqrt(L*C); fc=wn/(2*pi); xita=1/(2*wn*R*C); z=tf('z',Ta); Gpd=c2d(Gpc,Ta); F=Gpd.num; D=Gpd.den; n=F{:,:}; d=D{:,:}; b1=n(2); b2=n(3); a1=d(2); a2=d(3); % Parametros del controlador OSAP modificado p1=a1; p2=a2; q1=b1; q2=b2; p1m=(p2-a1*p1)/q1; p2m=-p1*a2/q1; q1m=1/q1; q2m=(p1*b1-q2)/q1; q3m=p1*b2/q1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cdc=4700e-6; Rdc=30; Rs=0.1; % Matrices de espacio de estado para simulación con carga Carga=-1/Rc/Cf; A=[ 0 -1/Lf 1/Cf Carga]; B=[1/Lf;0]; F=[0;-1/Cf]; C=[0 1]; D=0; [Gcc,Hcc]=c2d(A,B,Ta/p); [Gcc,Fcc]=c2d(A,F,Ta/p); % Matrices de espacio de estado para simulación sin carga Carga=0; A=[ 0 -1/Lf 1/Cf Carga]; B=[1/Lf;0]; F=[0;-1/Cf]; C=[0 1]; D=0; [Gc,Hc]=c2d(A,B,Ta/p); [Gc,Fc]=c2d(A,F,Ta/p); % Matriz de Normalização Tn=[1/Ibase 0 0 1/Vbase]; % Matrizes da Planta Normalizadas An=Tn*A*inv(Tn); Bn=Tn*B*Vbase; Fn=Tn*F*Ibase; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Discretización de la planta [G,H]=c2d(An,Bn,Ta); % Matrices muestreadas con el atraso de implementación digital Td=Ta; [eAtD,H1]=c2d(An,Bn,Td); [eATtD,H2]=c2d(An,Bn,Ta-Td); % Matriz H0 Ho=eATtD*inv(An)*(eAtD-eye(2))*Bn; % Matriz H1 H1=inv(An)*(eATtD-eye(2))*Bn; % MATRICES Gp y Hp en tiempo discreto que modelan el atraso de la implementación digital Gp=[G Ho;zeros(1,3)]; Hp=[H1;1]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Vectores de conmutacion y Matrices de decomposicion para la modulación del actuador V0=[0 0]'; V1=[1 0]'; V2=[0 1]'; V3=[1 1]'; M1= 1; M2=-1; M=[M1 M2]; SVector0=[0 0 0 0]; SVectorf=[1 0 0 1]; SVectors=[1 1 1 1]; % Inicialización de las variables utilizadas en la simulación S=1; T1=p/2; T2=p/2; ACTR0=SVector0(1,:); ACTRf=SVectorf(1,:); ACTRs=SVectors(1,:); kc=1:na*p; % Vector que indica la magnitud de las variables continuas con periodo de muestreo Ta/p ka=1:na; % Vector que indica la magnitud de las variables discretas con periodo de muestreo Ta % Vectores de la planta en tiempo discreto Xm=zeros(2,length(ka)); vc=zeros(1,length(ka)); uk=zeros(1,length(ka)); upi=zeros(1,length(ka)); uPD_F=zeros(1,length(ka)); uc=zeros(1,length(ka)); e=zeros(1,length(ka)); y=zeros(1,length(ka)); r=zeros(1,length(ka)); td=zeros(1,length(ka)); iL=zeros(1,length(ka)); urp=zeros(1,length(ka)); Se=zeros(1,length(ka)); Id=zeros(1,length(ka)); Id1=zeros(1,length(ka)); uffc=zeros(1,length(ka)); dt=zeros(2,length(ka)); uosap_m=zeros(1,length(ka)); Vrms=zeros(1,length(ka)); vs=0; % Vectores de la planta en tiempo continuo x=zeros(2,length(kc)); Io=zeros(1,length(kc)); Ior=zeros(1,length(kc)); upwm=zeros(length(kc),1); PWM1=zeros(length(kc),1); PWM2=zeros(length(kc),1); io=zeros(length(kc),1); t=zeros(length(kc),1); portadora=zeros(length(kc),1); Vdc=zeros(length(kc),1); Vdc(1:100,1)=278*ones(100,1); Vo=zeros(1,length(kc)); Veficaz=zeros(1,length(kc)); Switp=zeros(length(kc),2); Timer1=0; flag=1; sinal_port=1; ComparReg1=p/2; ComparReg2=p/2; km=3; i=0; Tf=1/fr; Npd=Tf/Ta; Npc=Tf/(Ta/p); MA=0; CR=1; kr=0.15; Na=2; % Ganancias del PD predictivo zca = 0.08528; K1a = -0.73507; K2a = zca*K1a; tf=nc/fr; tsin=0:Ta:tf; ref=sin(2*pi*fr*tsin); % figure; plot(tsin,ref); % Inicia Barra de Progresión perc=0.01; passo=0.01; wb = waitbar(0,'Aguarde, Simulando Inversor Monofasico...'); % Se inicia simulación: ciclo for principal que simula la planta en tiempo continuo % con el actuador PWM for k=2:np t(k)=k*Ta/p; % Vector de tiempo continuo Timer1=Timer1 + sinal_port; portadora(k)=Timer1; % Muestreo y cálculo de la acción de control dentro del procesador digital if flag==1 flag=0; km=km+1; i=i+1; td(km)=km*Ta; % Vector de tiempo discreto Xm(:,km)=[x(1,k)/Ibase;x(2,k)/Vbase]; iL(km)=Xm(1,km); vc(km)=Xm(2,km); Vrms(km)=sqrt(1/Npd*sum(vc(1:Npd).^2)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Vector de estado xe(:,km)=[iL(km);vc(km);uk(km-1)]; % Referencia de la tensión de salida deseada % r(km)=(311/Vbase)*sin(2*pi*fr*km*Ta); % if k<0.5*np % r(km)=(75/Vbase)*sin(2*pi*fr*km*Ta); % end % Salida de la planta y(km)=vc(km); % Erro de Tensao r(km)=ref(i); r(km+1) = ref(i+1); e(km)=r(km) - y(km); % Calculo del Controlador Repetitivo con acción integral if km>N Se(km)=e(km)+Se(km-N); urp(km)=kr*Se(km-N+Na); end % Calculo del Controlador Plug-in Repetitivo (Tomizuka) % if km>N % urp(km)=kr*e(km - N + Na) + Qr*urp(km - N); % end % Controlador OSAP uosap_m(km)=q1m*r(km+1) + p1m*y(km-1) + p2m*y(km-2) + q2m*uosap_m(km-1) + q3m*uosap_m(km-2); % uk(km)=uosap_m(km) + urp(km); % Controlador PD predictivo más feedforward de la referencia uPD_F(km)=K1a*e(km-1) + K2a*e(km-2) + r(km); uk(km)=uPD_F(km); % + urp(km); if MA==1 uk(km)=r(km); end % MODULAÇÃO SPACE VECTOR % Identificação do setor em que o vetor uc se encontra if uk(km)>0 S=1; else S=2; end % CALCULO DOS TEMPOS T1 e T2 T1=p*[M((S-1)+1)]*uk(km)*(Vbase/Vcc); T2=p-T1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MAQUINA DE ESTADOS:IMPLEMENTA A MODULAÇAO SPACE VECTOR ACTR0=SVector0(S,:); ACTRf=SVectorf(S,:); ACTRs=SVectors(S,:); end if Timer1==p/2 sinal_port=-1; end if Timer1==0 sinal_port=1; flag=1; end if Timer1= 2.005 %0.085 Ior(k)=x(2,k)/Rc*2.0; end if t(k) >= 2.065 %0.145 Ior(k)=x(2,k)/Rc/10; end % Retificador não controlado com filtro capacitivo if abs(x(2,k))>=Vdc(k) Io(k)=sign(x(2,k))*(abs(x(2,k))-Vdc(k))/Rs; Vdc(k+1)= Vdc(k) + Ta/p/Cdc*(( abs(x(2,k)) - Vdc(k) )/Rs - Vdc(k)/Rdc ); else Io(k)=0; Vdc(k+1)= Vdc(k) + Ta/p/Cdc*(-Vdc(k)/Rdc ); end % if k<=(np/2)+(ms*p/16) % Io(k)=0; % end if CR==1 Io(k)=Ior(k); end % Calculo del valor eficaz de la tensión de salida Vo(k)=x(2,k); % Veficaz(k)=sqrt(1/Npc*sum(Vo(1:Npc).^2)); x(:,k+1)= Gc*x(:,k) + Hc*upwm(k) + Fc*Io(k); % Atualiza a barra de progressão if fix((round(np))*perc)==k waitbar(perc) perc=perc+passo; end end close(wb) figure plot(t,upwm) axis([3.88 max(t) -Vcc-10 Vcc+10]) figure plot(td(1:km),e(1:km)) legend('error') axis([0 max(t) -1.0 1.0]) figure plot(td,uk) legend('Acción de Control') axis([3.88 max(td) -1.1 1.1]) figure plot(td,vc,'r') hold plot(td,Vrms,'k') plot(td,r(1:km),'m') legend('vo','RMS','ref') axis([3.88 max(td) -1.1 1.1]) figure plot(t,x(2,1:k),'r') hold plot(t,5*Io) grid legend('vo','io') axis([1.96 2.12 -400 400]) THD_tension_salida