clear all close all % Definicion de los parametros de simulacion fm=5000; % FRECUENCIA DE MUESTREO Tm=1/fm; % PERIODO DE MUESTREO f=50; % FRECUENCIA EN Hz DE LA TENSION SENOIDAL DE SALIDA T1=1/f; % PERIODO DE LA FUNDAMENTAL DE LA TENSION DE SALIDA Tf=20*T1; % TIEMPO TOTAL DE SIMULACION na=Tf/Tm; % NUMERO DE TOTAL DE MUESTRAS EN LA SIMULACION p=1; % NUMERO DE PUNTOS POR CADA PERIODO DE MUESTREO (PASO DE SIMULACION) np=na*p; % NUMERO TOTAL DE PUNTOS DE LA SIMULACION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PARAMETROS DE LA PLANTA. Lf=250e-6; % Indutancia del filtro de salida en mH Cf=125e-6; % Capacitancia del filtro de salida en uF Rc=30; % Resistencia de carga en Ohms % Valores base para normalización de la planta a valores en p.u. Vbase=311; % VALOR DE LA FUENTE DE TENSION CC DE ENTRADA EN VOLTIOS Ibase=Vbase/Rc; % VALOR BASE DE CORRIENTE PARA NORMALIZACION (VALOR CC NOMINAL) EN AMPERES w=2*pi*f; % Frecuencia angular de la tension senoidal de salida en rad/seg T=Tm; Vref=100; % Definición de las matrices de estado del sistema % x1=vo, x2=iL A=[-1/Rc/Cf 1/Cf; -1/Lf 0]; B=[0;1/Lf]; F=[-1/Cf;0]; C=[1 0]; D=0; % Determinación de la controlabilidad del sistema rango_planta=rank(ctrb(A,B)); % rango debe ser igual al orden del sistema % Matriz de Normalizacion Tn=[1/Vbase 0 0 1/Ibase]; % Matrizes de la Planta Normalizadas An=Tn*A*inv(Tn); Bn=Tn*B*Vbase; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Discretizacion de la planta para proyecto del controlador por realimentación de estados [Gn,Hn]=c2d(An,Bn,Tm); % Determinación de la controlabilidad del sistema rango_ps=rank(ctrb(Gn,Hn)); % rango debe ser igual al orden del sistema: planta + servo % Especificaciones de proyecto Mp=1/100; % Sobrepaso maximo Mp=exp(-pi.sigma/wd) ts=1e-3; % Tiempo de establecimiento 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 s1=-sigma+1i*wd; % Polos dominantes de lazo cerrado s2=-sigma-1i*wd; % complejos conjugados en el dominio continuo snd=-3*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); % znd=0.0001; % Polo no dominante en el dominio discreto % Definición de las matrices de la planta con el servo % G(circunfleja), H(circunfleja) Gc=[ Gn Hn; zeros(1,2) 0]; Hc=[0;0;1]; Cc=[C 0]; Dc=0; % Proyecto usando ubicación de polos polos=[z1 z2 znd]; Kplace=place(Gc,Hc,polos); Kssp=(Kplace+[zeros(1,2) 1])*inv([Gn-eye(2) Hn;C*Gn C*Hn]); K2=Kssp(1:2); % Ganancias para realimentacion de estados K1=Kssp(3); % Ganancia del servo % Conjunto de ganancias de Schuster y Cia. % K2=[2.7367 -3.7139]; K1=-2.2690; % Determinación de la controlabilidad del sistema rango_comp=rank(ctrb(Gc,Hc)); % rango debe ser igual al orden del sistema: planta + servo % Verifica si los autovalores se corresponden con los deseados autovalores=eig(Gc-Hc*Kssp); % 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 % Determinación de la controlabilidad del sistema rango_planta=rank(ctrb(A,B)); % rango debe ser igual al orden del sistema % Discretizacion de las matrices para simulación de la planta en tiempo continuo % con periodo de muestreo=T/p. Planta sin normalizar [Gs,Hs]=c2d(A,B,T/p); [Gs,Fs]=c2d(A,F,T/p); % Discretizacion de las matrices para simulación de la planta en tiempo continuo % con periodo de muestreo=T/p. Planta sin normalizar Rc=10; A=[-1/Rc/Cf 1/Cf; -1/Lf 0]; B=[0;1/Lf]; F=[-1/Cf;0]; C=[1 0]; D=0; [Gs1,Hs1]=c2d(A,B,T/p); [Gs1,Fs1]=c2d(A,F,T/p); % Variables de la planta discreta xd=zeros(2,length(ka)); xm=zeros(2,length(ka)); vo=zeros(1,length(ka)); iL=zeros(1,length(ka)); y=zeros(1,length(ka)); r=zeros(1,length(ka)); rampa=zeros(1,length(ka)); rd=zeros(1,length(ka)); e=zeros(1,length(ka)); v=zeros(1,length(ka)); t=zeros(1,length(ka)); u=zeros(1,length(ka)); uc=zeros(1,length(ka)); Wc=zeros(1,length(ka)); VO=zeros(1,length(ka)); IO=zeros(1,length(ka)); IC=zeros(1,length(ka)); IL=zeros(1,length(ka)); Flag_var=1; % 1 = Es para referencia sinusoidal % 0 = Es para variar referencia en escalón y verificar desempeńo Flag_carga=1; % 1 = Para para variación de carga for k=2:na t(k)=k*T; rampa(k)=t(k); % r(k)=311; r(k)=(311/Vbase)*sin(2*pi*f*t(k)); % Perfil de variación de referencia para verificar % especificaciones. Arranque en rampa. if k < na/(na/1000) && Flag_var == 1 rd(k) = (1/0.2)*rampa(k)*r(k); else rd(k)=(311/Vbase)*sin(2*pi*f*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(k)=Vref/Vbase; end if k > 0.5*na && Flag_var == 0 rd(k)=311/Vbase; end xm(:,k)=[xd(1,k)/Vbase;xd(2,k)/Ibase]; vo(k)=xm(1,k); iL(k)=xm(2,k); y(k)=vo(k); e(k)=rd(k)-y(k); v(k)=v(k-1) + e(k); u(k)=-K2*xm(:,k) + K1*v(k); uc(k)=Vbase*u(k); Gss=Gs; Hss=Hs; if k > 0.7*na+0.25*(T1/Tm) && Flag_carga == 1 Gss=Gs1; Hss=Hs1; end xd(:,k+1)=Gss*xd(:,k) + Hss*uc(k); VO(k)=xd(1,k); IC(k)=Cf*(VO(k)-VO(k-1))/T; IL(k)=xd(2,k); IO(k)=IL(k)-IC(k); end % figure; plot(t,rd) figure; plot(t,y,'r',t,rd,'k'); grid; legend('tension de salida'); figure; plot(t,e,'r'); grid; legend('error de tension'); axis([0 t(length(t)) -1 1]) figure; plot(t,v,'k'); grid; legend('servo') figure; plot(t,u,'m'); grid; legend('accion de control') figure; plot(t,VO); grid; legend('tension de salida') figure; plot(t,IL); grid; legend('corriente de el inductor') figure; plot(t,IC); grid; legend('corriente de el capacitor') figure; plot(t,VO,t,IO); grid; legend('tensión de salida','corriente en la carga')