% Control Digital y no Lineal. Departamento de Electronica - FI - UNAM % Responsable: Ing. Fernando Botterón. % Este programa simula en el dominio de tiempo discreto la aplicación del % Principio del Modelo Interno. En este caso se desea que la salida de un % convertidor CC-CA monofasico full-bridge con modulación PWM, cuya planta es modelada % como un sistema de segundo orden (Filtro pasa bajas LC), % siga asintoticamente una referencia senoidal de una frecuencia y amplitud % determinadas. Se presentan simulaciones con el modelo interno que usa el generador % de señales periódicas y como compensadores estabilizador, se utiliza o un PD predictivo o % una realimentación de los estados de la planta, pudiendo seleccionarse uno u otro % según se desee. El inversor monofasico puede ser simulado tanto con carga lineal (resistiva pura) % como con carga no lineal (rectificador a diodos con filtro capacitivo). % Finalmente, un script .m externo calcula la THD de la tensión sintetizada por este convertidor. % DEFINICION DE PARAMETROS DE SIMULACION clear close all format long f1=50; % FRECUENCIA DESEADA DE LA TENSION DE SALIDA EN Hz T1=1/f1; w=round(2*pi*f1); % FRECUENCIA ANGULAR DESEADA DE LA TENSION DE SALIDA EN rad/seg mf=200; % NUMERO DE PULSOS PWM POR PERIODO fsw=mf*f1; % FRECUENCIA DE CONMUTACION Tsw=1/fsw; % PERIODO DE CONMUTACION fm=1*fsw; % FRECUENCIA DE MUESTREO T=1/fm; % PERIODO DE MUESTREO fmi=fm; % FRECUENCIA DE MUESTREO DEL MODELO INTERNO Tmi=1/fmi; % PERIODO DE MUESTREO DEL MODELO INTERNO ms=round(fm/f1); % NUMERO DE MUESTRAS POR PERIODO DE LA TENSION DE SALIDA N=round(fmi/f1); % NUMERO DE MUESTRAS DEL MODELO INTERNO EN UN PERIODO DE LA TENSION DE SALIDA, % Y TAMBIEN EL ORDEN DEL CONTROLADOR POR MODELO INTERNO SIN INCLUIR UN POSIBLE FILTRO Q(z) nc=20; % NUMERO DE PERIODOS DE SIMULACION na=ms; % NUMERO DE TOTAL DE MUESTRAS POR SIMULACION npp=mf*nc; % NUMERO TOTAL DE PULSOS PWM Tclock=25e-09; % PERIODO DE CLOCK DEL TIMER INTERNO DEL DSP 25e-09 = 40,00 MHz TPWM=round(Tsw/Tclock); % VALOR EN CUENTAS PARA EL PERIODO DE CONMUTACION: Tsw/Tclock TPER=TPWM/2; % VALOR EN CUENTAS A CARGAR EN EL REGISTRO DE PERIODO PARA EL MODULO PWM: PULSO CENTRADO Tp=TPER/2; % VALOR EN CUENTAS DEL PERIODO PARA EL CALCULO DEL CICLO UTIL: TPWM/4 p=TPWM; % NUMERO DE PUNTOS POR PERIODO DE MUESTREO PARA DETERMINAR EL PASO DE SIMULACION np=npp*p; % NUMERO TOTAL DE PUNTOS DE LA SIMULACION % Valores base para la normalización de la planta a valores en p.u. para % facilitar la representación dentro del procesador digital Pbase=5000; % VALOR BASE DE POTENCIA EN VOLT-AMPERES Vbase=311; % VALOR BASE DE TENSION PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA TENSION DE SALIDA) EN VOLTS Ibase=32; % VALOR BASE DE CORRIENTE PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA CORRIENTE DE CARGA) EN AMPERES Vcc=400; % VALOR DE LA FUENTE DE TENSION CC DE ENTRADA EN VOLTS Vpico=311; % Valor de pico deseado de la tensión de salida. Puede elegirse otro valor si se desea una tensión eficaz diferente Rc=(Vbase/sqrt(2))/(Ibase/sqrt(2)); % Resistencia de carga para la potencia nominal en Ohms %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DE LA PLANTA. Filtro LC y Carga Resistiva Pura Lf=500e-6; % Indutancia en uH Cf=120e-6; % Capacitancia en uF wo=1/sqrt(Lf*Cf); % Frecuencia angular de resonancia del filtro LC en rad/seg fo=wo/2/pi; % Frecuencia de resonancia del filtro LC en Hz %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DEL RECTIFICADOR NO CONTROLADO UTILIZADO COMO CARGA NO LINEAL Cdc=470e-6; % Capacitancia del filtro en uF Rdc=15; % Resistencia de carga en Ohms Rs=0.5; % Resistencia serie para suavizar la di/dt de carga del capacitor Cdc % Matrices de estado de la planta % Vector de estados: x=[iL;vc]. iL=corriente a traves del inductor, vc=tension en bornes del capacitor A=[ 0 -1/Lf 1/Cf 0 ]; B=[1/Lf;0]; F=[0;-1/Cf]; C=[0 1]; D=0; Mtc=ctrb(A,B); estados_controlables=rank(Mtc); % Matriz de Normalizacion Tn=[1/Ibase 0 0 1/Vbase]; % Matrices de la Planta Normalizadas An=Tn*A*inv(Tn); Bn=Tn*B*Vbase; Fn=Tn*F*Ibase; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Discretizacion de la planta normalizada para proyecto del controlador por realimentación de estados [G,H]=c2d(An,Bn,T); % Obtención de las matrices de la planta discreta en espacio de estados, considerando el atraso de transporte % de la implementacion digital en tiempo real Td=T; % Se considera un tiempo de atraso igual al periodo de discretizacion T [eAtD,H1]=c2d(An,Bn,Td); [eATtD,H2]=c2d(An,Bn,T-Td); % Matriz H0 Ho=eATtD*inv(An)*(eAtD-eye(2))*Bn; % Matriz H1 H1=inv(An)*(eATtD-eye(2))*Bn; % MATRICES Gp y Hp DE LA PLANTA DISCRETA QUE TIENEN EN CUENTA EL ATRASO DE TRANSPORTE Gp=[G Ho;zeros(1,3)]; Hp=[H1;1]; Cp=[0 1 0]; Mtd=ctrb(Gp,Hp); estados_controlables=rank(Mtd); % Determinación de la controlabilidad del sistema rango_ps=rank(ctrb(Gp,Hp)); % rango debe ser igual al orden del sistema: planta + atraso digital d_lqr = 0; % Especificaciones de proyecto Mp=1/100; % Sobrepaso maximo Mp=exp(-pi.sigma/wd) ts=1e-3; % Tiempo de asentamiento ts=4/sigma sigma=4/ts; % Parte real de los polos dominantes de lazo cerrado wd=-pi*sigma/log(Mp); % Parte imaginaria de los polos dominantes de lazo cerrado Kdom1=3; Kdom2=3.001; s1=-sigma+1i*wd; % Polos dominantes de lazo cerrado s2=-sigma-1i*wd; % complejos conjugados en el dominio continuo snd1=-Kdom1*sigma; % Polo no dominante en el dominio continuo snd2=-Kdom2*sigma; % Polo no dominante en el dominio continuo z1=exp(s1*T); % Polos dominantes de lazo cerrado z2=exp(s2*T); % complejos conjugados en el dominio discreto znd1=exp(snd1*T); znd2=exp(snd2*T); % Polo no dominante en el dominio discreto % Definición de las matrices de la planta con el servo % G(circunfleja), H(circunfleja) Gc=[ Gp Hp; zeros(1,3) 0]; Hc=[0;0;0;1]; Cc=[Cp 0]; Dc=0; % Determinación de la controlabilidad del sistema rango_comp=rank(ctrb(Gc,Hc)); % rango debe ser igual al orden del sistema: planta + servo % Proyecto usando ubicación de polos polos=[z1 z2 znd1 znd2]; Kplace=place(Gc,Hc,polos); Kssp=(Kplace + [zeros(1,3) 1])*inv([Gp-eye(3) Hp;Cp*Gp Cp*Hp]); K2=Kssp(1:3); % Ganancias para realimentacion de estados K1=Kssp(4); % Ganancia del servo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Proyecto del controlador por realimentación de estados: Tecnica DLQR if d_lqr == 1 Q=[10 0 0 0 0 10 0 0 0 0 1 0 0 0 0 1]; R=1; Klqr=dlqr(Gc,Hc,Q,R); Kssp=(Klqr + [zeros(1,3) 1])*inv([Gp-eye(3) Hp;Cp*Gp Cp*Hp]); K2=Kssp(1:3); % Ganancias para realimentacion de estados K1=Kssp(4); % Ganancia del servo end %% Aqui se verifica el diseño a través de simulación del desempeño para entrada en escalon Gs =[ Gp-Hp*K2 Hp*K1 (Cp*Hp*K2 - Cp*Gp) 1 - Cp*Hp*K1]; Hs = Hc; sistema_pca = ss(Gs,Hs,Cc,Dc,T); tca = 0:T:0.02; figure; step(Vbase*sistema_pca,tca) % [y,t,x]=step(sistema_pca,tca); figure; pzmap(sistema_pca) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Discretizacion de las matrices para simulación de la planta en tiempo continuo con periodo de muestreo=T/p Tdp=T/TPWM; % 1- Discretizacion utilizando un muestreador-retenedor de orden cero (ZOH): comando "c2d" [Gsim,Hsim]=c2d(A,B,Tdp); [Gcc,Fsim]=c2d(A,F,Tdp); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % INICIALIZACION DE LAS VARIABLES UTILIZADAS EN LA SIMULACION kc=1:np; % Vector con total de puntos de simulación de las variables en tiempo continuo % con periodo de muestreo T/p ka=1:na; % Vector con total de puntos de simulación de las variables en tiempo discreto % con periodo de muestreo T % Variables de la planta discreta xm=zeros(2,length(ka)); xp=zeros(3,length(ka)); vc=zeros(1,length(ka)); u=zeros(1,length(ka)); un=zeros(1,length(ka)); e=zeros(1,length(ka)); y=zeros(1,length(ka)); r=zeros(1,length(ka)); rd=zeros(1,length(ka)); rampa=zeros(1,length(ka)); td=zeros(1,length(ka)); iL=zeros(1,length(ka)); umi=zeros(1,length(ka)); v=zeros(1,length(ka)); % Variables de la planta continua x=zeros(2,length(kc)); Io=zeros(1,length(kc)); ic=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); muestra=zeros(length(kc),1); CMP1=zeros(length(kc),1); CMP2=zeros(length(kc),1); Vdc=zeros(length(kc),1); Vdc(1:100,1)=278*ones(100,1); vo=zeros(1,length(kc)); % Valores iniciales de los contadores, comparadores e interrupciones Timer1=0; % Timer asociado al contador del módulo PWM: up-down flag=1; % Flag que simula interrupción del módulo PWM sign_counter=1; % signo del contador para generar la portadora up-down ComparReg1=TPWM/2; % valores iniciales de los comparadores de salida del modulo PWM ComparReg2=TPWM/2; % Inicializacion del indice de avance de la simulacion discreta, debido a % que Matlab no acepta indices cero ni negativos. km=1; LA=0; % flag de operacion en lazo abierto: 0 opera a lazo cerrado CR=0; % flag de operacion con carga resistiva: 0 opera con carga no lineal Flag_var=1; % 1 = Es para variar referencia y verificar desempeño % 0 = Es para referencia en escalón A=Vpico/Vbase; % Amplitud de la referencia % Factor de normalizacion para ancho de pulso: Define la profundidad de modulacion fn=(Vbase/Vcc); % Inicia Barra de Progresión perc=0.01; passo=0.01; wb = waitbar(0,'Aguarde, Simulando Inversor Monofasico...'); % Inicio de la Simulación for k=1:np % k es el indice de avance de la simulacion en tiempo continuo t(k)=k*Tsw/TPWM; % t(k) es el vector de tiempo continuo % Muestreo y calculo de la accion de control: % Aqui se produce la interrupcion del periferico del DSC o uC. Por % ejemplo del modulo PWM o del modulo ADC if flag==1 flag=0; km=km+1; % km es el indice de avance de la simulacion en tiempo discreto td(km)=(km-1)*T; % t(k) es el vector de tiempo discreto rampa(km)=t(k); % Referencia de la tension que se desea sintetizar a la salida del inversor. r(km)=A*sin(2*pi*f1*t(k)); % Perfil de variación de referencia para verificar % especificaciones. Arranque en rampa. if km < 0.2*nc*na && Flag_var == 1 rd(km) = (1/(0.2*nc*T1))*rampa(km)*r(km); else rd(km)=A*sin(2*pi*f1*t(k)); end % Fin perfil de variación de referencia en rampa% % Muestreo de la tension de salida y de la corriente en el inductor. % Aqui se almacenan en memoria RAM las variables de estado, que % seran utilizadas para el calculo del controlador por % realimentacion de estados. xm(:,km)=[x(1,k)/Ibase;x(2,k)/Vbase]; iL(km)=xm(1,km); vc(km)=xm(2,km); % El signo está negativo porque la polaridad de las señales PWM está invertida % Vector de estado de la planta con atraso de transporte % Aqui se define el vector de estado de la planta modelada en el espacio de estado que tiene en % cuenta el atraso Td de la conversion de las variables mas el calculo de la accion de control. xp(:,km)=[iL(km);vc(km);u(km-1)]; % Error de Tension e(km)=rd(km) - vc(km); v(km)=v(km-1) + e(km); % Accion de Control aplicada a la planta: Modelo Interno mas realimentacion de estados u(km)=K1*v(km) - K2*xp(:,km); if u(km) > 1 u(km) = 1; v(km) = (u(km) + K2*xp(:,km))/K1; end if u(km) < -1 u(km) = -1; v(km) = (u(km) + K2*xp(:,km))/K1; end if LA==1 % Operación en lazo abierto u(km)=rd(km); end % Calcula el ancho de pulso: Porcentaje de Tsw un(km)=u(km)*fn; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Define contador up-down para pulso PWM centrado if Timer1==TPWM/2 sign_counter=-1; flag=1; end if Timer1<0 sign_counter=1; % flag=1; end % Actualizacion de los registros comparadores con el valor del ancho de pulso % para comparacion con el contador up-down if Timer1==TPWM/2 % ComparReg1 = Tp + un(km)*Tp; ComparReg1 = TPER/2 + un(km)*(TPER/2); % ComparReg2 = Tp - un(km)*Tp; ComparReg2 = TPER/2 - un(km)*(TPER/2); end % Maquina de estado que define el patron PWM if ComparReg1 > Timer1 PWM1(k)=1; else PWM1(k)=0; end if ComparReg2 > Timer1 PWM2(k)=1; else PWM2(k)=0; end % Tension PWM aplicada a la planta upwm(k)=(PWM1(k)-PWM2(k))*Vcc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Corriente de carga lineal resistiva Ior(k)=x(2,k)/Rc/100; % Variacion de carga resistiva % Esto nos permite efectuar una variacion en escalon de la carga lineal % para verificar el rechazo a disturbio del sistema controlado % if k>=(0.41*np)+(ms/4) % Ior(k)=x(2,k)/Rc; % end % % if k>=(0.66*np)+(ms/4) % Ior(k)=x(2,k)/Rc/10; % end if t(k) >= 0.166 Ior(k)=x(2,k)/Rc; end if t(k) >= 0.266 Ior(k)=x(2,k)/Rc/100; end % A continuacion, se simula la carga no lineal: % Rectificador monofasico con diodos con filtro capacitivo y carga resistiva. 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) + T/p/Cdc*(( abs(x(2,k)) - Vdc(k) )/Rs - Vdc(k)/Rdc ); else Io(k)=0; Vdc(k+1)= Vdc(k) + T/p/Cdc*(-Vdc(k)/Rdc ); end % Si CR = 1 entonces el disturbio de carga es lineal: Carga resistiva. % En caso contrario, Io(k) será la corriente drenada por el rectificador. if CR==1 Io(k)=Ior(k); end vo(k)=x(2,k); %Tension de salida usada para el calculo de la THD ic(k)=x(1,k) - Io(k); % Corriente en el capacitor del filtro de salida % Calcula la ecuación a diferencias de la planta en tiempo continuo % obteniendo los estados del sistema dinamico (corriente y tension) % para la accion de control actual, upwm(k) y el disturbio de carga actual, Io(k) x(:,k+1)= Gsim*x(:,k) + Hsim*upwm(k) + Fsim*Io(k); % Actualiza la barra de progresion if fix((round(np))*perc)==k waitbar(perc) perc=perc+passo; end % Incremento del contador para comparación y generacion de PWM Timer1=Timer1 + sign_counter; portadora(k)=Timer1; muestra(k)=flag; CMP1(k) = ComparReg1; CMP2(k) = ComparReg2; end close(wb) % Cierra la barra de progresion % Graficos de las variables discretas figure plot(td(1:km),e(1:km)) legend('error') axis([0 max(td) -1 1]) figure plot(td(1:km),v(1:km)) legend('servo') % axis([0 max(td) -1 1]) figure plot(td(1:km),u(1:km)) legend('acción de control') % axis([0 max(td) -1 1]) figure plot(td(1:km),vc(1:km),'r') hold plot(td(1:km),rd(1:km),'k') grid legend('vc','referencia') % axis([0.18 max(td) -1.1 1.1]) % Graficos de las variables continuas figure plot(t(1:k),upwm,t(1:k),x(2,1:k),'r') % axis([max(t)-0.04 max(t) -380 380]) figure plot(t,x(2,1:k),'r') hold plot(t,1*Io) grid legend('vc','io'); % axis([0.18 0.2 -350 350]) figure plot(t,x(1,1:k),'r') hold plot(t,1*Io); plot(t,ic,'k') grid legend('iL','io','ic'); % axis([max(t)-0.04 max(t) -120 120]) % figure % plot(t,CMP1) % hold on % plot(t,CMP2,'r') % plot(td,un*Tp - Tp-1000,'k') % legend('CMP1','CMP2','un'); grid % % axis([max(t)-0.04 max(t) -2*Tp 2*Tp]) % % figure % subplot(2,1,1); plot(t(1:k),CMP1,'b',t,CMP2,'r') % title('Comparadores CMP1(k) - CMP2(k)'); % % axis([max(t)-0.03 max(t) 0 2*Tp]); % subplot(2,1,2); plot(td(1:km),un,'k'); % title('Acción de control Normalizada'); % % axis([max(t)-0.03 max(t) -1 1]); % % figure % plot(t,CMP1,'b',t,CMP2,'r',t,PWM1*Tp+3*Tp+2000,t,PWM2*Tp+2*Tp+500,t,portadora,'r',t,muestra*0.5*Tp/2,'k',t,(PWM1*Tp - PWM2*Tp) - Tp-2000,'r'); % % axis([max(t)-0.04 max(t) -4*Tp 5*Tp]); % legend('CMP1','CMP2','PWM1','PWM2','Timer','muestras','upwm'); % % figure % subplot(2,1,1); plot(t,CMP1,'b',t,CMP2,'r',t,portadora,'m',t,muestra*0.2*TPER,'k') % % axis([0.03 0.051 -200 2*Tp+200]); % title('Comparadores, Contador Up-Down y Muestreo'); % subplot(2,1,2); plot(t,upwm); % % axis([0.03 0.051 -400 400]); % title('Tension PWM'); % Calcula la THD y grafica el espectro de la tension de salida THD_tension_salida_HB