%=========UFSC - Federal University of Santa Catarina========
% Author: Tiago Davi Curi Busarello
% Date: August 2019
% Project: Highly Accurate Digital Current Controllers for Single-Phase 
% LCL-Filtered Grid-Connected Inverters
% The paper is open-source and it available on 
% https://www.mdpi.com/2673-4826/1/1/2
%
%
% Obs: This script has the design for 4 controller. You should enable and
% disabale the parts for each controller

clear all
% close all
format long
clc 

s = tf('s');
opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';

%---------------------------System Parameters----------------------------
Urms = 127; %Line-to-Line RMS voltage
Pn = 1000; %Rated Power
f0 = 60; %Grid frequency
w0 = 2*pi*f0; %angular grid frequency
fsw = 10e3; %Switching frequency
fs = fsw; %Sampling frquency
Ts = 1/fs; %Sampling period
Ta = Ts;
fa = fs;
Cdc = 3e-3; %Capacitance of the DC-link
Hi = 1/10;
VDC = 220;

%---------------------Controller Parameters------------------------------
fc = fs/8;  %desired cut-off frequency
MFd = 60;  %desired phase margin
ts = 1/fsw;
Vtri = 2;  %Amplitude of the carrier considering double update of the PWM (1/(2*Cpk))
Ginv = 2*VDC/Vtri; %Inverter gain

%-----------------------LCL paramteres-----------------------------------
Lc = 2.28e-3;
Lg = 9.9e-4;
Cf = 1.64e-6;
Rc = 0.01;
Rg = 0.01;

fres = sqrt((Lc+Lg)/(Lc*Lg*Cf))/(2*pi); %Resonant frequency
Rd = 1/(2*pi*fres*Cf); %Damping resistor

%Checking constrains
f_lower = 10*f0;
f_higher = fsw/2;

%-------------------------Transfer functions-----------------------------
HLCL = (s*Cf*Rd+1)/(Lc*Lg*Cf*s^3 + Cf*(Lc+Lg)*Rd*s^2 + (Lc+Lg)*s);

bode(HLCL,opts)
title('Frequency response of the designed LCL filter ');
% hold
% bode(LCLapprox,opts)
grid

td = 1/(1.5*fs);
Delay = -(s-(2/td))/(s+(2/td));

%------------------Design Type 3 controller with and without considering
%the PWM delay ---------------------------------------------------------

FTMAu = Ginv*HLCL*Hi*Delay; %Transfer function considering the delay
% FTMAu = Ginv*HLCL*Hi;   %Transfer function not considering the delay

figure
bode(FTMAu,opts)
title('Uncontrolled Open Loop Transfer Function');
% xlim([0.001 1e5]);
grid

inst = evalfr(FTMAu, 2*pi*fc*j);
fase_fc = rad2deg(phase(inst));

ganho_fc = mag2db(abs(inst));
ganho_real = 10^(abs(ganho_fc/20));

alpha = MFd - fase_fc - 90;
K2 = (tand(alpha/4 + 45))^2;
R1 = 1000;
C2 = 1/(2*pi*fc*ganho_real*R1);
C1 = C2*(K2-1);
R2 = sqrt(K2)/(2*pi*fc*C1);
R3 = R1/(K2-1);
C3 = 1/(2*pi*fc*R3*sqrt(K2));

Cv = (s^2*R2*C1*C3*(R1+R3) + s*(R2*C1+R1*C3+R3*C3)+1) /(s^3*R1*R2*R3*C1*C2*C3 + s^2*(R1*R3*C3*(C1+C2)+R1*R2*C1*C2)+s*R1*(C1+C2));

FTMAc = FTMAu*Cv;
figure
bode(FTMAu,opts);
title('Bode for uncontrolled transfer function');
hold 
bode(FTMAc,opts)
title('Bode for controlled TF with Type 3 controller');
bode(Cv,opts)
title('Bode of designed Type 3 controller');
grid

FTMF = FTMAc/(FTMAc + 1);
figure 
bode(FTMF,opts)
title('Bode for closed-loop with type 3 controller');
grid

%Discretização do tipo 3
D1 = R2*C1*C3*(R1+R3);
D2 = R2*C1 + R1*C3 + R3*C3;
D3 = R1*R2*R3*C1*C2*C3;
D4 = R1*R3*C1*(C1+C2) + R1*R2*C1*C2;
D5 = R1*(C1+C2);

Dcte = D3+D4*ts+D5*ts^2;

b0 = (ts*D1 + D2*ts^2 + ts^3)/Dcte;
b1 = -(2*ts*D1 + D2*ts^2)/Dcte;
b2 = D1*ts/Dcte;
b3 = 0;

a0 = 1;
a1 = -(3*D3+2*D4*ts+D5*ts^2)/Dcte;
a2 = (3*D3+D4*ts)/Dcte;
a3 = -D3/Dcte;

b = [b0 b1 b2 b3];
a = [a0 a1 a2 a3];

H = tf(b,a,ts);

%------------------Design Type 2 controller without considering PWM delay

% FTMAu = Ginv*HLCL*Hi;   %Transfer function not considering the delay
% 
% figure
% bode(FTMAu,opts)
% title('Uncontrolled Open Loop Transfer Function');
% xlim([0.001 1e5]);
% grid
% 
% inst = evalfr(FTMAu, 2*pi*fc*j);
% fase_fc = rad2deg(phase(inst));
% 
% ganho_fc = mag2db(abs(inst));
% ganho_real = 10^(abs(ganho_fc/20));
% 
% alpha = MFd - fase_fc - 90;
% K2 = tand(alpha/2 + 45);
% 
% Prod = 1/(2*pi*fc*ganho_real*K2);
% C2 = 100e-9;
% R1 = Prod/C2;
% C1 = C2*((K2^2)-1);
% R2 = K2/(2*pi*fc*C1);
% 
% Ctipo2 = (1+s*R2*C1)/(s*R1*(C1+C2+s*R2*C1*C2));
% 
% figure
% bode(Ctipo2,opts)
% title('Bode for designed Type 2 controller');
% xlim([0.001 1e5]);
% grid
% 
% FTMAc = FTMAu*Ctipo2;
% figure
% bode(FTMAc,opts)
% xlim([0.001 1e5]);
% title('Bode for Controlled Transfer Function with Type 2 controller');
% grid
% 
% FTMF = FTMAc/(FTMAc + 1);
% figure 
% bode(FTMF,opts)
% title('Bode for closed-loop with type 2 controller');
% grid
% 
% %Discretization of type 2 controller
% 
% beta = 4*R1*R2*C1*C2 + 2*ts*R1*(C1+C2);
% 
% b0 = (ts^2 + 2*ts*C1*R2)/beta;
% b1 = (2*ts^2)/beta;
% b2 = (ts^2 - 2*ts*C1*R2)/beta;
% 
% a0 = 1;
% a1 = (-8*R1*R2*C1*C2)/beta;
% a2 = (4*R1*R2*C1*C2 - 2*ts*R1*(C1+C2))/beta;
% 
% b = [b0 b1 b2];
% a = [a0 a1 a2];
% 
% H = tf(b,a,ts);

%------------------Design PR controller without considering PWM delay

% FTMAu = Ginv*HLCL*Hi;   %Transfer function not considering the delay
% 
% figure
% bode(FTMAu,opts)
% title('Uncontrolled Open Loop Transfer Function');
% xlim([0.001 1e5]);
% grid
% 
% %Desired Parameters of the Proportional-Resonant Controller
% fr = 60; %Desired resonant frequency [Hz]
% wr = 2*pi*fr; %Desired resonant angular frequency [rad/s]
% qsi = 0.95; %Damping factor
% Bs = 1.5; %Desired resonant frequency bandwidth [Hz]
% Br = 2*pi*Bs; %Desired resonant angular frequency bandwidth [rad/s]
% 
% %Design of the proportional and resonant gains
% Linv = Lc+Lg;
% Rinv = Rc + Rg;
% alpha = 2*qsi + 1;
% Kp = (alpha*sqrt(alpha)*wr*Linv-Rinv)/((VDC/2)*Hi); %Desgined proportional gain
% Ki = (wr^2)*Linv*((alpha^2)-1)/(VDC*Hi); %Deesigned resonant gain
% 
% %Design of the resonant filter
% kr = 1; %Gain of the resonant filter - not including the Ki
% Cte = ((kr*(Br^2)*0.5)/sqrt((wr^2)-0.25*(Br^2)))*(exp(-0.5*Br*Ta))*sin(Ta*sqrt((wr^2)-0.25*(Br^2))); %Defined constant for simplicity
% 
% b0 = kr*Br*Ta;
% b1 = (-kr*Br*(exp(-0.5*Br*Ta))*cos(Ta*sqrt((wr^2)-0.25*(Br^2)))-Cte)*Ta;
% b2 = 0;
% 
% a0 = 1;
% a1 = -2*(exp(-0.5*Br*Ta))*cos(Ta*sqrt((wr^2)-0.25*(Br^2)));
% a2 = exp(-Br*Ta);
% 
% b = [b0,b1,b2];
% a = [a0,a1,a2];
% 
% TFz = tf(b,a,ts); %Defining the transfer function in z-domina
% z = (2+Ta*s)/(2-Ta*s); %Bilinear transformation for mapping z-domain in w-domain
% 
% TFw = (b0 + b1*z^-1 + b2*z^-2)/(a0 + a1*z^-1 + a2*z^-2); %Transfer function in w-domain
% TFPRw = TFw*Ki+Kp;
% % 
% % figure
% % bode(TFw,opts)
% % title('Frequency response in w-domain of the resonant part');
% % grid
% % 
% figure
% bode(TFPRw,opts)
% title('Frequency response in w-domain of the PR controller');
% grid
% 

% ------------Design of the active damping coeffcient and verifying -------

[Gm] = margin(FTMAc);
Ka_min = Gm*2*pi*fc*Lc/VDC;
Ka_max = 4*Vtri*fsw*Lc/VDC;

Ka = (Ka_min+Ka_max)/2;  %Active damping coefficient

T = Ginv*Hi/(Lc*Lg*Cf*s^3 + Cf*Lg*Ka*Ginv*s^2 + (Lc+Lg)*s);

% Verify from cut-off frequency and lower if they are coincident. If yes,
% the controller can be designed considering the transfer function of
% passive damping. 

figure
bode(FTMAu,opts)
hold
bode(T,opts)
title('Comparing the transfer functions considering and not considering the damping factor. The passive and active damping are coincident for regions for f<fc');
grid

%To adjust phase offset: Right click --> properties --> options --> adjust
%phase offset. 