% 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 f1=50; % FRECUENCIA DESEADA DE LA TENSION DE SALIDA EN Hz w=round(2*pi*f1); % FRECUENCIA ANGULAR DESEADA DE LA TENSION DE SALIDA EN rad/seg mf=; % 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=; % VALOR BASE DE POTENCIA EN VOLT-AMPERES Vbase=; % VALOR BASE DE TENSION PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA TENSION DE SALIDA) EN VOLTS Ibase=; % VALOR BASE DE CORRIENTE PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA CORRIENTE DE CARGA) EN AMPERES Vcc=; % VALOR DE LA FUENTE DE TENSION CC DE ENTRADA EN VOLTS Vpico=; % 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=; % Indutancia en uH Cf=; % 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=1000e-6; % Capacitancia del filtro en uF Rdc=10; % 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=corriente a traves del inductor, vc=tension en bornes del capacitor % Matriz de Normalizacion Tn=[1/Xbase 0 0 1/Xbase]; % Matrices de la Planta Normalizadas An=Tn*A*inv(Tn); Bn=Tn*B*Xbase; Fn=Tn*F*Xbase; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Proyecto del controlador por realimentación de estados: Reubicación de polos % Determinación de la controlabilidad y de la observabilidad del sistema % discreto rango_mc=; % rango debe ser igual al orden de la planta % Especificaciones de proyecto Mp=; % Sobrepaso maximo Mp=exp(-pi.sigma/wd) ts=; % Tiempo de establecimiento ts=4/sigma s1=-sigma+1i*wd; % Polos dominantes de lazo cerrado s2=-sigma-1i*wd; % complejos conjugados en el dominio continuo snd=-x*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 znd=exp(snd*T); % Definición de las matrices de la planta con el servo % G(circunfleja), H(circunfleja) % Proyecto usando ubicación de polos polos=; Kplace=; Kssp=; K2=Kssp(1:2); % Ganancias para realimentacion de estados K1=Kssp(3); % Ganancia del servo % Determinación de la controlabilidad de la planta + el servo controlador rango_mcs=rank(ctrb(Gc,Hc)); % rango debe ser igual al orden del sistema a lazo cerrado: planta + servo %% Aqui se verifica el diseño a través de simulación del desempeño para entrada en escalon Gs=[ G-H*K2 H*K1 (C*H*K2 - C*G) 1 - C*H*K1]; Hs=Hc; sistema_pca=ss(Gs,Hs,Cc,Dc,T); tca=0:T:6; figure; step(sistema_pca,tca) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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" [Gc,Hc]=c2d(A,B,Tdp); [Gcc,Fc]=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)); Se=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=3; % Flags para definir las diferentes formas de operacion y diferentes tipos de carga LA=1; % flag de operacion en lazo abierto: 0 opera a lazo cerrado CR=1; % 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 Flag_carga=0; % 1 = Para para variación de carga Qmi=1; % Filtro Q escalar para el controlador por modelo interno A=Vpico/Vbase; % Amplitud de la referencia % Factor de normalizacion para ancho de pulso: Define la profundidad de modulacion fn=(Vbase/Vcc); Vref=155; % 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-3)*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 % % Perfil de variación de referencia en escalón if Flag_var == 0 rd(km)=Vref/Vbase; end if k > 0.5*na && Flag_var == 0 rd(km)=Vpico/Vbase; end % Muestreo de las variables de estado xm(:,km)=[x(1,k)/Xbase;x(2,k)/Xbase]; x1(km)=xm(1,km); x2(km)=xm(2,km); % reemplazar x1 y x2 por las variables físicas correspondientes % 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)=[x1(km);x2(km);u(km-1)]; y(km)=Cp*xp(:,km); %% Ecuaciones de cálculo dentro del procesador % Ecuacion del error e(km)=; % Ecuacion del servo controlador % Ecuacion de la accion de control u(km)=; if LA==1 % Operación en lazo abierto u(km)=r(km); end % Calcula el ancho de pulso: Porcentaje de Tsw un(km)=u(km)*fn; % Acción de control normalizada 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; % 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 Flag_carga == 1 Ior(k)=x(2,k)/Rc/10; if k>=(0.6*np)+(ms*p/4) Ior(k)=x(2,k)/Rc; end if k>=0.85*np Ior(k)=x(2,k)/Rc/10; end 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)= Gc*x(:,k) + Hc*upwm(k) + Fc*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),vc(1:km),'r') hold plot(td(1:km),rd(1:km),'k') grid legend('vc','referencia') axis([0 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