Contents

% 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

Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0     273    4.776062e+03    1.562e+00    4.027e+02
    1     546    2.535763e+03    4.706e-01    2.522e+02    3.722e+00
    2     819    1.699587e+03    6.374e-02    1.027e+02    2.806e+00
    3    1092    1.484663e+03    0.000e+00    5.127e+01    1.239e+00
    4    1365    1.240195e+03    0.000e+00    9.293e+00    1.627e+00
    5    1638    8.967204e+02    0.000e+00    1.029e+01    3.230e+00
    6    1911    7.978252e+02    0.000e+00    6.674e+00    4.179e-01
    7    2184    6.935036e+02    7.823e-03    7.258e+00    7.075e-01
    8    2457    6.446517e+02    8.329e-03    8.158e+00    6.637e-01
    9    2730    6.202406e+02    4.479e-03    8.296e+00    7.400e-01
   10    3003    6.068349e+02    9.207e-04    8.776e+00    9.380e-01
   11    3276    6.039324e+02    1.098e-03    8.595e+00    7.592e-01
   12    3549    6.033423e+02    3.354e-04    8.869e+00    5.119e-01
   13    3823    6.021314e+02    0.000e+00    9.763e+00    7.981e-01
   14    4096    5.995387e+02    3.409e-03    8.725e+00    9.214e-01
   15    4369    5.970058e+02    8.611e-03    1.702e+01    1.115e+00
   16    4642    5.970251e+02    4.980e-03    5.501e+00    2.553e-01
   17    4916    5.975300e+02    1.096e-04    2.962e+00    9.095e-01
   18    5190    5.978709e+02    0.000e+00    2.968e+00    4.565e-01
   19    5464    5.978475e+02    0.000e+00    1.584e+00    2.107e-01
   20    5737    5.917980e+02    3.175e-03    1.589e+00    6.434e-01
   21    6010    5.893784e+02    2.967e-03    1.608e+00    6.480e-01
   22    6285    5.892908e+02    2.220e-03    1.605e+00    1.165e-01
   23    6558    5.890096e+02    3.490e-04    1.567e-01    3.550e-01
   24    6831    5.882545e+02    1.061e-03    1.553e-01    3.486e-01
   25    7104    5.881029e+02    3.640e-04    1.478e-01    1.726e-01
   26    7379    5.880920e+02    2.592e-04    1.155e-01    2.340e-02
   27    7652    5.880569e+02    0.000e+00    8.526e-03    9.545e-02
   28    7925    5.879277e+02    9.784e-05    4.844e-02    1.217e-01
   29    8198    5.878866e+02    5.658e-05    3.062e-02    6.611e-02
   30    8471    5.878693e+02    2.112e-05    3.034e-02    3.817e-02

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
   31    8744    5.878645e+02    6.556e-06    1.572e-02    1.592e-02
   32    9017    5.878624e+02    1.550e-05    1.584e-02    2.007e-02
   33    9290    5.878625e+02    8.295e-06    4.149e-03    6.356e-03
   34    9564    5.878626e+02    0.000e+00    4.086e-03    1.602e-02
   35    9838    5.878626e+02    0.000e+00    4.092e-03    7.267e-03
   36   10112    5.878626e+02    0.000e+00    2.489e-03    3.774e-03
   37   10385    5.878626e+02    0.000e+00    5.552e-04    1.569e-03
   38   10658    5.878626e+02    0.000e+00    5.510e-04    7.242e-04
   39   10931    5.878626e+02    0.000e+00    1.729e-04    9.931e-04
   40   11204    5.878623e+02    0.000e+00    7.089e-04    8.038e-04
   41   11478    5.878623e+02    0.000e+00    7.087e-04    6.819e-04
   42   11752    5.878623e+02    0.000e+00    2.729e-04    1.995e-04
   43   12026    5.878623e+02    0.000e+00    1.247e-04    9.839e-05
   44   12299    5.878623e+02    0.000e+00    6.871e-05    5.240e-05

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.




exitflag =

     1

Elapsed time is 1035.826782 seconds.

Published with MATLAB® R2018a