% Team Members: Sherin Adel Helal , Omniyah Gul M Khan
% Topic: "A New Protection Scheme Considering Fault Ride Through Requirements for Transmission Level Interconnected Wind Parks"
% Abstract: In this project, we attempt to re-simulate the results obtained by Saleh et al. using their proposal of a communication
% based TCV DOCR tripping characteristics. An IEEE 24-bus transmission system with  DFIG based wind turbine data was obtained and 
% was used to develop the protection coordination optimization model using the proposed Dual TCV approach to obtain the optimal 
% relay settings. Results were compared with the ITC method used in the literature to compare the FRT violations and DOCRs 
% tripping times. The ITC method was formulated as a linear programming problem and MATLABs simplex technique was used to obtain 
% the optimal time dial settings (TDS) of each relay.  The Dual-TCV method was formulated as a non-linear programming problem and 
% was solved using the sequential quadratic programming algorithm to obtain the optimal time dial settings (TDS) of each relay as
% well as the third relay setting (K) that controls the rate of change in the tripping time with respect to change in fault voltage magnitude. 

%% In this m file we simulate the main Dual TCV Method

clc;clear all; warning off;format long;
tic
A = [];
b = [];
Aeq = [];
beq = [];
load('Results24bus_TCV.mat')
lb = [0.05*ones(1,68);-inf*ones(1,68);0.05*ones(1,68);-inf*ones(1,68)];

x0 = [0.05*ones(1,68);1*ones(1,68);0.05*ones(1,68);1*ones(1,68)];

ub = [1*ones(1,68);inf*ones(1,68);1*ones(1,68);inf*ones(1,68)];

options = optimset('Display','iter','MaxIter',100000,'MaxFunEvals',10000000);
%fmincon- NLP solver.
%fmincon(objective function,initial values of opt. variable, "Ax =< B where
%A & B relate to our inequality constraints","AeqX =< Beq where
%A & B relate to our equality constraints",lower bound on opt. variable X, Upper bound on opt. Variable X, NL function C(x) targeting inequality constraints
%options= as set using optimset to decide on max number of iterations and
%function evaluations).
[x,fval,exitflag,output,lambda,grad,hessian]=fmincon(@myfun_24busp,x0,A,b,Aeq,beq,lb,ub,@const_fun_24busp,options);

%Usage of optimal solutions 
TDS1 = x(1,:);
K1 = x(2,:);
TDS2 = x(3,:);
K2 = x(4,:);

load('Ip.mat');
Ip=Ipp;
load('Fault_analysis_24bus_8dfig.mat')

first_fl=1;
last_fl=29;
time_sum=0;


for i=1:4 %i: counter denoting fault types (1:SLG, 2:DLG, 3:LL,4: 3phase SC)
    for j=1:3 %j: fault location- nearend, midpoint, farend of line.
        
        clear ISC V  t ISCb ISCc Vb Vc
        
        if i==1 || i==4
            ISC=Fault_analysis.(['ISCa' num2str(i) num2str(j)]);
            V=Fault_analysis.(['Vfa' num2str(i) num2str(j)]);
        elseif i==2 || i==3
            
            ISCb=Fault_analysis.(['ISCb' num2str(i) num2str(j)]);
            ISCc=Fault_analysis.(['ISCc' num2str(i) num2str(j)]);
            Vb=Fault_analysis.(['Vfb' num2str(i) num2str(j)]);
            Vc=Fault_analysis.(['Vfc' num2str(i) num2str(j)]);
            ISC=max(ISCb,ISCc);
            V=min(Vb,Vc);
        end
        flx=1;
        for fl=first_fl:last_fl
            
            MAP=MAPPING(fl,first_fl);    %This function gives us the list of primary and backup relays for each fault location
            for k=1:2       %k represents the two primary relays on each fault location
                indp(k)=MAP(k,1);   %the primary relay number
                indb1(k)=MAP(k,2);  %the backup for the kth primary relay
                indb2(k)=MAP(k,3);  %the second backup for the kth primary relay
                indb3(k)=MAP(k,4);  %the third backup for the kth primary relay
                indb4(k)=MAP(k,5);  %the fourth backup for the kth primary relay
                              
                
                if indp(k)~=0   %Primary relay tripping time
                    
                    t(indp(k),fl)=((1/(exp(1-V(indp(k),fl))))^K1(indp(k)))*0.14*TDS1(indp(k))/(-1 + (ISC(indp(k),fl)/Ip(indp(k)))^.02 );
                    time_sum=time_sum+ t(indp(k),fl);
                    if k==1 %primary relay 1 at fault location (fl) 
                        All_time.(['t' num2str(i) num2str(j)])(2*flx-1,1)=t(indp(k),fl);
                    else %primary relay 2 at fault location (fl) -> looking from opp. direction to Primary relay 1
                        All_time.(['t' num2str(i) num2str(j)])(2*flx,1)=t(indp(k),fl);
                    end
                end
                
                
                if indb1(k)~=0 %Backup 1 relay tripping time
                    
                    t(indb1(k),fl)=((1/(exp(1-V(indb1(k),fl))))^K2(indb1(k)))*0.14*TDS2(indb1(k))/(-1 + (ISC(indb1(k),fl)/Ip(indb1(k)))^.02 );
                    time_sum=time_sum+ t(indb1(k),fl);
                    if k==1 %backup1 relay 1 at fault location (fl) 
                        All_time.(['t' num2str(i) num2str(j)])(2*flx-1,2)=t(indb1(k),fl);
                    else %backup1 relay 2 at fault location (fl) -> looking from opp. direction to backup1 relay 1
                        All_time.(['t' num2str(i) num2str(j)])(2*flx,2)=t(indb1(k),fl);
                    end
                end
                
                
                
                if indb2(k)~=00 %Backup2 relay tripping time

                    t(indb2(k),fl)=((1/(exp(1-V(indb2(k),fl))))^K2(indb2(k)))*0.14*TDS2(indb2(k))/(-1 + (ISC(indb2(k),fl)/Ip(indb2(k)))^.02 );
                    time_sum=time_sum+ t(indb2(k),fl);
                    if k==1 %backup2 relay 1 at fault location (fl)
                        All_time.(['t' num2str(i) num2str(j)])(2*flx-1,3)=t(indb2(k),fl);
                    else  %backup2 relay 2 at fault location (fl) -> looking from opp. direction to backup2 relay 1
                        All_time.(['t' num2str(i) num2str(j)])(2*flx,3)=t(indb2(k),fl);
                    end
                end
                
                
                if indb3(k)~=0%Backup3 relay tripping time
                    
                    
                    t(indb3(k),fl)=((1/(exp(1-V(indb3(k),fl))))^K2(indb3(k)))*0.14*TDS2(indb3(k))/(-1 + (ISC(indb3(k),fl)/Ip(indb3(k)))^.02 );
                    time_sum=time_sum+ t(indb3(k),fl);
                    if k==1 %backup3 relay 1 at fault location (fl)
                        All_time.(['t' num2str(i) num2str(j)])(2*flx-1,4)=t(indb3(k),fl);
                    else  %backup3 relay 2 at fault location (fl) -> looking from opp. direction to backup3 relay 1
                        All_time.(['t' num2str(i) num2str(j)])(2*flx,4)=t(indb3(k),fl);
                    end
                end
                
               
                if indb4(k)~=0 %Backup4 relay tripping time        
                    t(indb4(k),fl)=((1/(exp(1-V(indb4(k),fl))))^K2(indb4(k)))*0.14*TDS2(indb4(k))/(-1 + (ISC(indb4(k),fl)/Ip(indb4(k)))^.02 );
                    time_sum=time_sum+ t(indb4(k),fl);
                    if k==1 %backup4 relay 1 at fault location (fl)
                        All_time.(['t' num2str(i) num2str(j)])(2*flx-1,5)=t(indb4(k),fl);
                    else  %backup4 relay 2 at fault location (fl) -> looking from opp. direction to backup4 relay 1
                        All_time.(['t' num2str(i) num2str(j)])(2*flx,5)=t(indb4(k),fl);
                    end
                end
                
                
            end
            flx=flx+1;
        end
        
        Op_time.(['t' num2str(i) num2str(j)])= t;
        
        
        
    end
end
save('Results24bus_TCV.mat','All_time','Op_time', 'TDS1','K1', 'TDS2','K2','time_sum');
exitflag
toc