% 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 half-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=100; % NUMERO DE PULSOS PWM POR PERIODO (25 para 1250 Hz). Cambiar a 100 para 5kHz fsw=mf*f1; % FRECUENCIA DE CONMUTACION Tsw=1/fsw; % PERIODO DE CONMUTACION fm=2*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 FUNDAMENTAL 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=30; % NUMERO DE PERIODOS DE SIMULACION na=ms; % NUMERO DE TOTAL DE MUESTRAS POR PERIODO DE 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 (DSP TMS320LF2407) 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=10000; % 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=64.28; % VALOR BASE DE CORRIENTE PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA CORRIENTE DE CARGA) EN AMPERES Vcc=2*350; % VALOR DE LA FUENTE DE TENSION CC DE ENTRADA EN VOLTS: (311 para indice de modulación = 1) % Aumentar el valor a 350V para operación a lazo cerrado 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=1500e-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=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;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; % 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]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Proyecto del controlador por realimentación de estados: Tecnica DLQR Q=[1 0 0 0 1 0 0 0 1]; R=1; Kre=dlqr(Gp,Hp,Q,R); % Ganancias del compensador PD predictivo obtenidas con el script "PD_predictivo_solve.m" K1=-0.02; K2=0.01; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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: variables internas que maneja el DSP 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)); 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=p/2; % valores iniciales de los comparadores de salida del modulo PWM ComparReg2=p/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, diferentes % tecnicas de control y diferentes tipos de carga LA=0; % flag de operacion en lazo abierto: 0 opera a lazo cerrado PD=1; % flag de operacion con PD predictivo RE=0; % flag de operacion con realimentacion de estados CR=1; % flag de operacion con carga resistiva: 0 opera con carga no lineal 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 limitando % el pulso mínimo y el pulso máximo para lograr que el inversor opere siempre en la región lineal fn=(Vbase/Vcc); % Definicion de la constante del controlador por modelo interno y del % avance para compensacion de la estabilidad en las altas frecuencias % Para operacion con realimentacion de estados if RE==1 kmi=0.2; d=3; end % Para operacion con proporcional-derivativo predictivo if PD==1 kmi=0.15; d=3; end % 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 % 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); % 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-1);vc(km-1);u(km-1)]; % Referencia de la tension que se desea sintetizar a la salida del inversor. r(km)=A*sin(2*pi*f1*t(k)); % Error de Tension e(km)=r(km) - vc(km); % Calculo del Controlador por Modelo Interno: Utiliza la cadena de % integradores con realimentacion positiva con el filtro Q(z) if km>N umi(km)=kmi*e(km - N + d) + Qmi*umi(km - N); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if RE==1 % Accion de Control aplicada a la planta: Modelo Interno mas realimentacion de estados ure(km)=umi(km) - Kre*xp(:,km); u(km)=ure(km); end if PD==1 % Accion de Control aplicada a la planta: Modelo Interno mas accion PD predictiva % No es necesaria la suma de la accion feedforward [o sea:r(km)] upd(km)= K1*e(km-1) + K2*e(km-2); % + r(km); u(km)=1*umi(km) + upd(km); end if LA==1 % Operación en lazo abierto u(km)=r(km); end % Normaliza la acción de control según el valor base y la tensión % CC de entrada al convertidor un(km)=u(km)*fn; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Define contador up-down para pulso PWM centrado if Timer1<0 sign_counter=1; flag=1; end if Timer1>TPWM/2 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 = un(km)*2*Tp + Tp; % Puede usarse indistintamente una ecuación ComparReg1 = un(km)*TPER + TPER/2; % u otra end if Timer1==0 % ComparReg1 = un(km)*2*Tp + Tp; ComparReg1 = un(km)*TPER + TPER/2; end % Maquina de estado que define el patron PWM if ComparReg1 > Timer1 PWM1(k)=1/2; else PWM1(k)=-1/2; end % Tension PWM aplicada a la planta upwm(k)=PWM1(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 k>=(0.2*np)+(ms*p/4) % Ior(k)=x(2,k)/Rc*2; % end % % if k>=0.5*np % Ior(k)=x(2,k)/Rc/10; % 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 de la simulación 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; % Genera la portadora triangular para graficarla muestra(k)=flag; % Guarda las muestras para graficarlas CMP1(k) = ComparReg1; % Guarda los valores del comparador para graficar 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),r(1:km),'k') grid legend('vc','referencia') axis([max(t)-0.04 max(t) -1.1 1.1]) % Graficos de las variables continuas figure plot(t,upwm,t,x(2,1:k),'r') axis([max(t)-0.04 max(t) -Vcc/2 Vcc/2]) figure plot(t,x(2,1:k),'r') hold plot(t,1*Io) grid legend('vc','io'); axis([max(t)-0.04 max(t) -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) -110 110]) figure plot(t,CMP1) hold on plot(td,un*2*Tp + Tp,'k') legend('CMP1','un'); grid axis([max(t)-0.02 max(t) 0 2*Tp]); figure subplot(2,1,1); plot(t(1:k),CMP1,'b') title('Comparador'); axis([max(t)-0.02 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.02 max(t) -0.5 0.5]); figure plot(t,CMP1,'m',t,portadora,'r',t,muestra*0.5*(p/2),'k',t,PWM1*Tp-Tp); axis([max(t)-0.02 max(t) -2*Tp 2*Tp+2000]); legend('CMP1','Contador','Muestreo','PWM'); figure subplot(2,1,1); plot(t,CMP1,'m',t,portadora,'r',t,muestra*0.25*TPER,'k') axis([max(t)-0.02 max(t) -500 2*Tp+500]); title('Comparador, Contador Up-Down y Muestreo'); subplot(2,1,2); plot(t,PWM1*2*Tp); axis([max(t)-0.02 max(t) -Tp-500 Tp+500]); title('Tension PWM'); % Calcula la THD y grafica el espectro de la tension de salida THD_tension_salida_HB