%=========UFSC - Federal University of Santa Catarina========
% Author: Tiago Davi Curi Busarello
% Date: June 2018
% Project: Design of a Digital Proportional-Resonant Current Controller
%
% Simulation of the paper Design Procedure for a Digital 
%   Proportional-Resonant Current Controller in a Grid Connected Inverter 
%   presented in 4th IEEE Southern Power Electronics Conference, SPEC 2018

close all
clear all
clc
format long

%Sampling Paramters
fa = 30000; %Sampling frequency in Hz
Ta = 1/fa;  %Samping period
wa = 2*pi*fa; %Sampling angular frequency

%Grid-Connected Inverter and System Parameters
Linv = 10e-3; %Output filter inductor [H]
Rinv = 0.5e-3; %Output filter resistance [ohm]
Vdc = 450; %Total DC-Link Voltage [V]
Vp = 180; %Peak voltage of the grid [V]
fg = 60; %Grid frequency [Hz]
wg = 2*pi*fg; %Angular grid frequency [rad/s]
Hi = 1/10; %Current sensor gain [A/A]
fsw = 30e3; %Switching frequency [Hz]
Pnom = 1000; %Inverter nominal power [W]

%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
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,fa); %Defining the transfer function in z-domina
s = tf('s');
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;

opts = bodeoptions('cstprefs');
opts.FreqUnits = 'Hz';
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

%Desired Parameters of the Proportional-Resonant Controller for the second
%resonant part - subscript 2 means second controller

fr2 = 300; %Desired resonant frequency [Hz]
wr2 = 2*pi*fr2; %Desired resonant angular frequency [rad/s]
qsi2 = 0.95; %Damping factor
Bs2 = 1.5; %Desired resonant frequency bandwidth [Hz]
Br2 = 2*pi*Bs2; %Desired resonant angular frequency bandwidth [rad/s]

%Design of the proportional and resonant gains
alpha2 = 2*qsi2 + 1;
Kp2 = (alpha2*sqrt(alpha2)*wr2*Linv-Rinv)/((Vdc/2)*Hi); %Desgined proportional gain
Ki2 = (wr2^2)*Linv*((alpha2^2)-1)/(Vdc*Hi); %Deesigned resonant gain

%Design of the resonant filter
kr2 = 1; %Gain of the resonant filter - not including the Ki
Cte2 = ((kr2*(Br2^2)*0.5)/sqrt((wr2^2)-0.25*(Br2^2)))*(exp(-0.5*Br2*Ta))*sin(Ta*sqrt((wr2^2)-0.25*(Br2^2))); %Defined constant for simplicity

b02 = kr2*Br2*Ta;
b12= (-kr2*Br2*(exp(-0.5*Br2*Ta))*cos(Ta*sqrt((wr2^2)-0.25*(Br2^2)))-Cte2)*Ta;
b22 = 0;

a02 = 1;
a12 = -2*(exp(-0.5*Br2*Ta))*cos(Ta*sqrt((wr2^2)-0.25*(Br2^2)));
a22 = exp(-Br2*Ta);

b2 = [b02,b12,b22];
a2 = [a02,a12,a22];

TFz2 = tf(b2,a2,fa); %Defining the transfer function in z-domina

TFw2 = (b02 + b12*z^-1 + b22*z^-2)/(a02 + a12*z^-1 + a22*z^-2); %Transfer function in w-domain
TFPRw2 = TFw2*Ki2+Kp2;

TF = TFw + TFw2;

figure
bode(TF,opts)
title('Frequency response in w-domain of the resonant parts combined');
grid

open_system('DigitalPR.slx')









