% 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 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 en conexión % paralela. Este modelo interno usa el generador de señales periódicas visto en clase. % Además, 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 all 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=80; % NUMERO DE PULSOS PWM POR PERIODO fsw=mf*f1; % FRECUENCIA DE CONMUTACION Tsw=1/fsw; % PERIODO DE CONMUTACION fm=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=nc*ms; % NUMERO DE TOTAL DE MUESTRAS POR SIMULACION p=512; % NUMERO DE PUNTOS POR CADA PERIODO DE MUESTREO (PASO DE SIMULACION) np=na*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 Vbase=311; % VALOR BASE DE TENSION PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA TENSION DE SALIDA) Ibase=50; % VALOR BASE DE CORRIENTE PARA NORMALIZACION (VALOR DE PICO MAXIMO DE LA CORRIENTE DE CARGA) Vcc=400; % VALOR DE LA FUENTE DE TENSION CC DE ENTRADA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DE LA PLANTA. Filtro LC y Carga Resistiva Pura Lf=400e-6; % Indutancia en uH Cf=700e-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 Rc=Vbase/Ibase; % Resistencia de carga para la potencia nominal en Ohms %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DEL RECTIFICADOR NO CONTROLADO UTILIZADO COMO CARGA NO LINEAL Cdc=1000e-6; % Capacitancia del filtro en uF Rdc=5; % 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; % Discretizacion de las matrices para simulación de la planta en tiempo continuo con periodo de muestreo=T/p Tdp=T/p; % 1- Discretizacion utilizando un muestreador-retenedor de orden cero (ZOH): comando "c2d" [Gc,Hc]=c2d(A,B,Tdp); [Gc,Fc]=c2d(A,F,Tdp); % 1- Discretizacion utilizando un muestreador-retenedor de orden cero (ZOH): % Gcd = transformada inversa de Laplace de la matriz [(sI - A)^-1] Gcd=[ cos(wo*Tdp) -wo*Cf*sin(wo*Tdp); wo*Lf*sin(wo*Tdp) cos(wo*Tdp)]; % Hcd=A^-1[Gcd - I]B Hcd=[wo*Cf*sin(wo*Tdp);(1-cos(wo*Tdp))]; % Hcd=A^-1[Gcd - I]F Fcd=[(1-cos(wo*Tdp));-wo*Lf*sin(wo*Tdp)]; % 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); z = tf('z',T); sysd = ss(G,H,C,D,T); % Transforma a de espacio de estado a FT [Num,Den] = tfdata(sysd,'v'); % en plano z, como num. y denom. Gpd = tf(Num,Den,T); % Función transferencia de la planta en z % Transfer function: % 0.09475 z + 0.09195 0.09475 + 0.09195 z^-1 b1 + b2 z^-1 % Gpd = ---------------------- = ------------------------ = ----------------- % z^2 - 1.728 z + 0.9146 z - 1.728 + 0.9146 z^-1 z + a1 + a2 z^-1 % Obtención del compensador % ------------------------------------------------------------------------ b1 = Num(2); b2 = Num(3); % Asigna coeficientes a1 = Den(2); a2 = Den(3); % Asigna coeficientes q1 = b1; q2 = b2; % Coeficientes del compensador. p1 = a1; p2 = a2; % Coeficientes del compensador. % 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.55; K2=0.01; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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)); uc=zeros(1,length(ka)); ure=zeros(1,length(ka)); upd=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)); uosap=zeros(1,length(ka)); Se=zeros(1,length(ka)); % Variables de la planta continua 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); muestra=zeros(length(kc),1); Vdc=zeros(length(kc),1); Vdc(1:100,1)=278*ones(100,1); vo=zeros(1,length(kc)); % Flags Timer1=0; flag=1; sinal_port=1; ComparReg1=p/2; 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 MA=0; % flag de operacion en lazo abierto PD=1; % flag de operacion con PD predictivo RE=0; % flag de operacion con realimentacion de estados CR=0; % flag de operacion con carga resistiva OS=0; % flag de operacion con OSAP perturba=1; % flag de variacion de carga Qmi=1; % Filtro Q escalar para el controlador por modelo interno A=1; % Amplitud de la referencia % 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.25; d=3; end % Para operacion con proporcional-derivativo predictivo if PD==1 kmi=0.2; d=3; end % Si fuera necesario, habria que agregar un kmi y un Na para operacion con % OSAP if OS==1 kmi=0.15; d=3; end % 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*T/p; % 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*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*km*T); % Error de Tension e(km)=r(km) - vc(km); % Calculo del Controlador por Modelo Interno: Utiliza la cadena de % integradores con realimentacion positiva + el filtro Qmi(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)= Kre*xp(:,km); u(km) = umi(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) = umi(km) + upd(km); end if OS==1 uosap(km)=(1/q1)*r(km) + (p1/q1)*vc(km) + (p2/q1)*vc(km-1) - (q2/q1)*uosap(km-1); u(km) = umi(km) + uosap(km); end if MA==1 % Operación en lazo abierto u(km)=r(km); end % Calcula el ancho de pulso: Porcentaje de Tsw uc(km)=u(km)*(p/4)*fn; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Define contador up-down para pulso PWM centrado if Timer1==p/2 sinal_port=-1; %flag=1; end if Timer1<0 sinal_port=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==0 ComparReg1= uc(km)+p/4; ComparReg2=-uc(km)+p/4; 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/10; % 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 perturba == 1 if k>=(0.2*np) Ior(k)=x(2,k)/Rc; end end % if k>=0.5*np % Ior(k)=x(2,k)/Rc; % 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 % 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 + sinal_port; portadora(k)=Timer1; muestra(k)=flag; 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.3 0.4 -1.1 1.1]) figure plot(td(1:km),vc(1:km),'r') hold plot(td(1:km),r(1:km),'k') legend('vc','referencia') axis([0.3 0.4 -1.1 1.1]) % Graficos de las variables continuas figure plot(t,upwm,t,x(2,1:k),'r') axis([0.3 0.321 -380 380]) if RE == 1 figure plot(td(1:km),u(1:km),td(1:km),ure(1:km),td(1:km),umi(1:km)) legend('Accion de control total','Accion de control RE','Accion de control MI') end if PD == 1 figure plot(td(1:km),u(1:km),td(1:km),upd(1:km),td(1:km),umi(1:km)) legend('Accion de control total','Accion de control PD','Accion de control MI') end if OS == 1 figure plot(td(1:km),u(1:km),td(1:km),uosap(1:km),td(1:km),umi(1:km)) legend('Accion de control total','Accion de control OSAP','Accion de control MI') end figure plot(t,x(2,1:k),'r') hold plot(t,2*Io) grid legend('vc','io'); axis([0.3 0.4 -350 350]) % figure % plot(t,-PWM1*0.5+0.55,t,-PWM2*0.5,t,(portadora/(p/2))+1.1,'r',t,muestra*0.5+1.1,'k'); axis([0.3 0.3008 -0.8 2.5]) % Calcula la THD y grafica el espectro de la tension de salida THD_tension_salida