%% Diseño de un control de tensión de salida de una fuente reductora de tensión
%% comúnmente conocida como convertidor buck. Se realiza el análisis en frecuencia
%% y en el tiempo de los efectos de la realimentación negativa en el control
%% de la tensión de salida y el rechazo a perturbaciones
clc
clear 
close all

s=tf('s');
P = bodeoptions; 
P.FreqUnits = 'Hz'; 
P.grid='on';

%% Definición de los parámetros de la planta 
%% Convertidor CC-CC reductor 
Po=2.5;         % Potencia de salida: 2.5 W
Vin=10;         % Tensión nominal de entrada en voltios
Von=2.5;        % Tensión nominal de salida en voltios
Io=Po/Von;      % Corriente nominal de carga
Ro=Po/Io^2;     % Resistencia de carga 
C=220e-6;       % Capacitancia del filtro de salida
L=75e-6;        % Inductancia del filtro de salida
rL=100e-3;      % Resistencia equivalente serie del inductor L
rC=70e-3;       % Resistencia equivalente serie del capacitor C
Vsw=2;          % Tensión de pico de la señal portadora diente de sierra en voltios
Ka=Vin/Vsw;     % Ganancia del actuador PWM
Fpwm=100e3;     % Frecuencia de conmutación de las llaves de potencia 100 kHz
Tpwm=1/Fpwm;    % Periodo de conmutación de la tensión Vpwm del actuador 10 us
Ks=1;           % Ganancia del sensor de tensión de salida
Vref=Von*Ks;    % Valor de referencia en voltios para operación a LC
Dn=Von/Vin;     % Ciclo útil nominal para buck ideal que define el punto de operación

%% Función de transferencia de la planta sin componentes de pérdidas
% Gp=Ka*Ro/(C*L*Ro*s^2 + L*s + Ro);
% zpk(Gp)

%% Función de transferencia de la planta: con RSE del inductor y del capacitor
% Gp=minreal(Ka*Ro*(rC*C*s + 1)/(s^2*C*L*(Ro + rC) + s*(L + rC*Ro*C + rL*C*(Ro + rC)) + Ro + rL));
% zpk(Gp)

%% Función de transferencia de la planta: con RSE del capacitor, rL = 0
Gp=minreal(Ka*Ro*(1 + rC*C*s)/(s^2*(C*L*Ro + C*L*rC) + (L + C*Ro*rC)*s + Ro));
zpk(Gp)

%% También puede escribirse la FT de la planta en base al factor Q
wz=1/(rC*C); fz=wz/2/pi; wo=1/sqrt(L*C); fo=wo/2/pi; Q=Ro*sqrt(C/L); H0=Vin/Vsw;
H=H0*(1 + s/wz)/((s/wo)^2 + s/(Q*wo) + 1);
% H=H0/((s/wo)^2 + s/(Q*wo) + 1);
zpk(H)

[Zp,Pp,Kplanta] = zpkdata(H,'v');

figure; bode(Gp,H,P); grid
legend('Rta de Gp','Rta de H(s)')

% figure; bode(H,P)
% legend('Rta de Gp')
%% Compensador
%% Ganancia o compensador proporcional
Gc=50;
%% PID: Cancela el cero de la planta
% Gc=(5*(s+30000)*(s+500))/(s*(s+6.494e04));

Gla=minreal(Gp*Gc*Ks);
Glc=minreal(feedback(Gla,Ks));

figure; bode(Gp,Gla,P)
figure; margin(Gla,P)
legend('Gla Sin Compensar','Gla Compensada')

figure; bode(Vref*Glc,P)  % Se multiplica por Vref para obtener el valor real en régimen
legend('Glc Compensada')

%% Respuesta de Gla y Gla/(1 + Gla)
figure; bode(Gla,Glc,P)
legend('Gla','Gla/(1 + Gla)')

%% Respuesta de 1/(1 + Gla)
S=minreal(1/(1 + Gla));
figure; bode(S,P)
legend('S = 1/(1 + Gla)')

figure; bode(S,Gla,P)
legend('1/(1 + Gla)','Gla')

%% Impedancia de salida de LA ideal sin pérdidas
Zoi= Ro*L*s/(C*L*Ro*s^2 + L*s + Ro);
figure; bode(Zoi,P)
legend('Zo LA sin perdidas')

%% Impedancia de salida de LA con ESR del inductor y del capacitor
Zoc=minreal((C*L*Ro*C*s^2 + (L*Ro + C*Ro*rC*rL)*s + Ro*rL)/((C*L*Ro + C*L*rC)*s^2 + (L + C*rC*rL + C*Ro*rC + C*Ro*rL)*s + Ro + rL));
figure; bode(Zoc,P)
legend('Zo LA con ESR del Inductor y del Capacitor')

%% Impedancia de salida de LA con ESR del inductor
ZoL=minreal((Ro*(s*L + rL))/(C*L*Ro*s^2 + (L + C*Ro*rL)*s + rL + Ro));
figure; bode(ZoL,P)
legend('Zo LA con ESR del Inductor')

%% Impedancia de salida de LA con ESR del capacitor
ZoC=minreal((L*Ro*s*(C*rC*s + 1))/((C*L*Ro + C*L*rC)*s^2 + (L + C*Ro*rC)*s + Ro));
figure; bode(ZoC,P)
legend('Zo LA con ESR del Capacitor')

%% Impedancia de lazo cerrado
Zlc=minreal(ZoC*S);
figure; bode(Zlc,P)
legend('Impedancia de salida de LC')

figure; bode(ZoC,Zlc,P)
legend('Impedancia de salida de LA','Impedancia de salida de LC')

%% Función de transferencia entre Vo y la tensión de entrada Ve a lazo abierto
Gpe=minreal(Vin*Ro*(1 + rC*C*s)/(s^2*(C*L*Ro + C*L*rC) + (L + C*Ro*rC)*s + Ro));
figure; bode(Gpe,P)
legend('FT Vo/Ve')

%% Función de transferencia entre Vo y la tensión de entrada Ve a lazo cerrado
Gvog=minreal(Gpe*S);
figure; bode(Gpe,Gvog,P)
legend('FT Vo/Ve','Relación de rechazo de entrada de LC')

paso_sim=1e-7; tfinal=4e-3;
t=0:paso_sim:tfinal-paso_sim;  % Vector de tiempo para graficar respuestas

close all

%% Respuesta a lazo abierto para el ciclo útil nominal
% figure; step(Dn*H);  % Permite observar la dinámica natural de la planta a LA

%% Respuesta a lazo abierto y a lazo cerrado sin compensador
Glasc=Gp*Ks;
figure; step(Glasc,t)
legend('Rta de LA SC: Vo sensor')

Glc_sc=feedback(Glasc,1);
figure; step(Vref*Glc_sc)
legend('Rta de LC SC: Vo sensor')

%% Respuesta a lazo cerrado compensada

Gu=feedback(Gc,Gp*Ks);
figure; step(Vref*Gu)
legend('Acción de control')

Ge=feedback(1,Gc*Gp*Ks);
figure; step(Vref*Ge)
legend('Error de tensión')

figure; step(Vref*Glc)
hold on
step(Vref*Glc_sc)
legend('Rta de LC: Vo sensor')

%% Las respuestas trazadas hasta aqui se corroboran tambien en el PSIM con los diagramas
%% de bloques y el circuito del convertidor con modulación PWM
% Error de posición:
% Kpos=dcgain(Glasc);
% essp=100/(1 + Kpos);

% Lugar de las raíces del sistema sin compesar
% figure; rlocus(Glasc);
% legend('LGR del sistema SC')
% close all

%% Diseño del controlador de adelanto por el lugar geométrico de las raíces
% Especificaciones de desempeño
% Mp=2/100;               % Sobrepaso maximo del 2%. Mp=exp(-pi.sigma/wd)
% ts=0.5e-3;              % Tiempo de establecimiento 0.5 mseg. ts=4.5/sigma
% sigma=4.5/ts;           % Parte real de los polos dominantes de lazo cerrado
% wd=-pi*sigma/log(Mp);   % Parte imaginaria de los polos dominantes de lazo cerrado
% 
% s1d=sigma+1i*wd;        % Polos dominantes de lazo cerrado
% s2d=sigma-1i*wd;        % complejos conjugados en el dominio continuo
% 
% pds=-sigma+1i*wd;       % Punto deseado
% Gdes=1/((s + s1d)*(s + s2d));
% 
% %% Se muestran el LGR incluyendo el polo al origen y especificaciones de desempeño
% %% para el control con un compensador integral
% % sisotool(Ka*Ks*Gp,(1/s),1,Vref)
% % figure; rlocus(Ka*Ks*Gp/s)
% %% Con un compensador integral no pueden cumplirse las especificaciones de desempeño
% %% dado que aunque se use una ganancia elevada, puede cumplirse el Mp pero no el ts.
% 
% % Polos y ceros de la FT de LA sin compensar con el polo al origen del integrador 
% 
% [Z,P,Kla] = zpkdata(Glasc/s,'v');
% p1=P(1);
% p2=abs(P(3));
% p3=abs(P(2));
% figure; pzplot(Glasc/s); hold on
% pzplot(Gdes); %axis([-4e4 0.5e4 -2.0e4 2.0e4])
% legend('Polos-ceros FT LA SC','Polo deseado s1d')
% 
% def=180-angle(evalfr(Glasc/s,pds))*180/(pi);
% angz=(180/pi)*angle(pds)/2 + def/(2);
% angp=(180/pi)*angle(pds)/2 - def/(2);
% 
% cero=wd/(tand(angz))+sigma;
% polo=wd/(tand(angp))+sigma;
% Kc=1/abs(evalfr(Glasc*(s+cero)/(s*(s+polo)),pds));
% 
% %% Compensador de adelanto 
% % Gc = Kc*(s + cero)/(s*(s + polo))
% Gc=Kc*(s + cero)/(s*(s + polo));
% 
% % A continuación se define la FT de LA compensada
% Gla=Gc*Glasc;
% figure; rlocus(Gla)
% legend('LGR resultante con Compensador de Adelanto e Integrador')
% % Para observar LGR con región deseada
% % sisotool(Ka*Gp,Gc,Ks,Vref)
% 
% % Se define la FT de LC compensada con compensador de adelanto de fase e integrador. 
% % En esta FT se obtiene la respuesta de la tensión de salida medida después del sensor.
% 
% Glcc=feedback(Gla,1);
% 
% tfinal=2e-3;
% t=0:paso_sim:tfinal-paso_sim;  % Vector de tiempo para graficar respuestas
% 
% % Obtención de las respuestas al escalón. Contrastar con las respuestas
% % obtenidas en buck_sincronico_LC_CC_Adelanto.psimsch
% % Se obtiene la respuesta al escalón de la tensión de salida del sensor
% figure; step(Vref*Glcc,t);
% legend('Vos - con Compensador de Adelanto e Integrador');
% 
% % Se obtiene la respuesta al escalón de la tensión de salida del convertidor
% figure; step(Vref*feedback(Gc*Ka*Gp,Ks),t);
% legend('Vo - con Compensador de Adelanto e Integrador');
% 
% %% Cálculo de la acción de control
% % Se define la FT entre la acción de control y la Vref
% Gr2u=feedback(Gc,Glasc);
% % Se obtiene la respuesta al escalón de la acción de control
% figure; step(Vref*Gr2u,t);
% legend('Acción de control para entrada en escalón')
% 
% % [y]=step(Vref*Glcc,t);
% % rstep=Vref*ones(length(t),1);
% % error=rstep-y;
% % figure; lsim(Gc,error,t);
% % legend('Acción de control para entrada en escalón')
% 
% %% Estabilidad del sistema compensado con PI
% figure; margin(Gla)
% legend('Margenes de Estabilidad del Buck')
% 
% figure; bode(Glcc)
% legend('Análsis de Desempeño en Frecuencia')
% 
% close all
% %% Se hace un reajuste del controlador para cumplir con el tiempo de asentamiento
% % sisotool(Ka*Gp,Gc,Ks,Vref)
% cero=9911; polo=23720; Kc=39615;
% Gcr=Kc*(s + cero)/(s*(s + polo));
% T=1/cero; alfa=1/T/polo;
% 
% % sisotool(Ka*Gp,Gcr,Ks,Vref)
% Glar=Gcr*Glasc;
% 
% figure; step(Vref*feedback(Gcr*Ka*Gp,Ks));
% legend('Vo - con Compensador de Adelanto e Integrador');
% 
% figure; margin(Glar)
% legend('Margenes de Estabilidad del Buck')
% 
% Glc=feedback(Gcr*Ka*Gp*Ks,1);
% P = bodeoptions; 
% P.FreqUnits = 'Hz'; 
% P.grid='on';
% 
% %% Respuesta de Gla y Gla/(1 + Gla)
% figure; bode(Glar,Glc,P)
% legend('Gla','Gla/(1 + Gla)')
% 
% %% Respuesta de 1/(1 + Gla)
% S=1/(1 + Glar);
% figure; bode(S,P)
% legend('S = 1/(1 + Gla)')
% 
% %% Respuesta de ZoLC
% figure; bode(Zo*S,P)
% 
% %% Cálculo de la acción de control
% % Se define la FT entre la acción de control y la Vref
% Gr2u=feedback(Gcr,Glasc);
% figure; step(Vref*Gr2u);
% legend('Acción de control para entrada en escalón')
% 
% %% Cálculo de los componentes del circuito electrónico
% % Compensador Adelanto: Ver presentación de clase unidad 4
% % T=R1*C1  alfa*T=R2*C2  Kc=R4*R2/(R3*R1)
% C1=150e-9;              % 150 nF
% C2=120e-9;              % 100 nF
% R1=T/C1;                % 1 kOHM    
% R2=alfa*T/C2;           % 360 OHM
% R3=100; Kc_alfa=Kc*alfa;
% R4=Kc_alfa*R3*R1/R2;    % 470 kOHM
% 
% close all
% %% Diseño de un controlador PID por cancelación polo-cero
% z1pid=abs(P(2)); z2pid=abs(P(3)); %z2pid=9685; 
% Mp=2/100;               % Sobrepaso maximo del 2%. Mp=exp(-pi.sigma/wd)
% ts=1.5e-3;              % Tiempo de establecimiento 0.5 mseg. ts=4.5/sigma
% sigma=4.5/ts;           % Parte real de los polos dominantes de lazo cerrado
% wd=-pi*sigma/log(Mp);   % Parte imaginaria de los polos dominantes de lazo cerrado
% 
% pds=-sigma+1i*wd;       % Punto deseado
% % Utilizo la condición de magnitud para calcular la ganancia del PID
% Kpid=1/abs(evalfr(Glasc*((s + z1pid)*(s + z2pid))/s,pds));
% % Kpid=97.8e-6;
% Gpid=Kpid*((s + z1pid)*(s + z2pid))/s;
% 
% % Para observar LGR con región deseada
% % sisotool(Ka*Gp,Gpid,Ks,Vref)
% 
% Glapid=Gpid*Glasc;
% figure; rlocus(Glapid)
% legend('LGR resultante con PID')
% 
% % A continuación se define la FT de LC compensada con PI
% % En esta FT se obtiene la respuesta de la tensión de salida
% % medida después del sensor.
% Glccpid=minreal(feedback(Glapid,1));
% 
% figure; margin(Glapid)
% 
% % Obtención de las respuestas al escalón. Contrastar con las respuestas
% % obtenidas en buck_sincronico_LC_CC_PI_AMPOP.psimsch
% figure; step(Vref*Glccpid,t);
% legend('Vos - Buck compensado con PID');
% 
% figure; step(Vref*feedback(Gpid*Ka*Gp,Ks),t);
% legend('Vo - Buck compensado con PID');
% 
% %% Cálculo de los componentes del circuito electrónico
% % PID con 1 solo AMP-OP: Ver presentación de clase unidad 4
% % Kc=(R1*C1)*(R2*C2)/(R1*C2)
% % tau_a=R1*C1; tau_b=R2*C2;
% format long
% tau_a=1/z1pid; tau_b=1/z2pid;
% % Seleccionamos C2
% C2=4.7e-6;
% % De la relacion de Kc obtenemos R1
% R1=tau_a*tau_b/(Kpid*C2);
% % De la relacion de tau_a obtenemos C1
% C1=tau_a/R1;
% % Y de la relacion de tau_b se obtiene R2
% R2=tau_b/C2;
% %% **************************************************************
% 
% %% Diseño del control ON-OFF con histéresis: Esto no va en el TP3
% % % La tensión de salida en regimen estable es
% % Vo = Von;
% % % La tensión a la salida del sensor en regimen estable es
% % Vs = Hs*Vo;
% % % A la vez, la tensión en la que se centra la ventana de histéresis es
% % Vcentro = Vs;
% % % La tensión de alimentación del comparador es
% % Vcc = 15;
% % % Para una ondulación del 2% se tiene que
% % DVo = 0.02*Vo;
% % % El voltaje de histéresis es un 50% de DVo
% % Vh = (0.5*DVo)*Hs;
% % % Los voltajes de nivel superior e inferior de basculación son:
% % VAS = Vcentro + Vh/2;
% % VAI = Vcentro - Vh/2;
% % % El factor n es:
% % n = (Vcc/Vh) - 1;
% % % La tensión de referencia es:
% % Vref = (((n + 1)/n)*Vcentro) - (Vcc/(2*n));
% % 
% % % Selecciono R2
% % R2 = 100;  % Ohms
% % R1 = n*R2;  % Ohms
% %% Adopto un valor comercial de 27k. 
% %% Esto hace que se obtenga un DVo = 0.205 V menor al nominal
% %% La ventana de histeresis resulta menor también con Vh = 0.055 V