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 comparison ITC Method

clc;clear all; warning off;format long;
tic
A = [];
b = [];
Aeq = [];
beq = [];

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

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

ub = [1*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_24bus,x0,A,b,Aeq,beq,lb,ub,@const_fun_24bus,options);
%Usage of optimal solution
TDS = x(1,:);

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   t ISCb ISCc

        if i==1 || i==4
            ISC=Fault_analysis.(['ISCa' 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)]);
            ISC=max(ISCb,ISCc);
        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)=0.14*TDS(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)=0.14*TDS(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)~=0 %Backup2 relay tripping time

                    t(indb2(k),fl)=0.14*TDS(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)=0.14*TDS(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 backup2 relay 1
                         All_time.(['t' num2str(i) num2str(j)])(2*flx,4)=t(indb3(k),fl);
                    end

                end

                if indb4(k)~=0

                    t(indb4(k),fl)=0.14*TDS(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_TC.mat','All_time','Op_time', 'TDS','k');
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      69    2.692857e+03    1.575e+00    2.699e+02
    1     138    2.101026e+03    1.019e+00    2.140e+02    8.871e-01
    2     207    2.034679e+03    9.266e-01    1.190e+02    1.012e-01
    3     276    2.008288e+03    8.516e-01    1.076e+02    4.665e-02
    4     345    1.942296e+03    6.102e-01    8.714e+01    1.331e-01
    5     414    1.940115e+03    5.226e-01    6.378e+01    3.327e-02
    6     483    1.961537e+03    3.163e-01    4.387e+01    8.740e-02
    7     552    2.016413e+03    2.697e-01    3.516e+01    1.348e-01
    8     621    2.077468e+03    2.356e-01    2.480e+01    1.253e-01
    9     690    2.208528e+03    2.078e-01    2.042e+01    2.515e-01
   10     759    2.259293e+03    1.757e-01    1.567e+01    9.531e-02
   11     828    2.269994e+03    1.708e-01    1.356e+01    1.984e-02
   12     897    2.357329e+03    1.301e-01    1.274e+01    1.694e-01
   13     966    2.433846e+03    9.342e-02    9.948e+00    1.491e-01
   14    1035    2.528511e+03    5.108e-02    9.027e+00    1.803e-01
   15    1104    2.547184e+03    4.065e-02    4.507e+00    4.162e-02
   16    1173    2.574275e+03    2.784e-02    1.688e+00    5.403e-02
   17    1242    2.594308e+03    1.896e-02    4.211e-02    3.822e-02
   18    1311    2.636085e+03    9.757e-04    2.570e-02    7.838e-02
   19    1380    2.636383e+03    5.585e-04    6.676e-03    1.445e-03
   20    1449    2.636786e+03    2.118e-04    1.513e-03    1.184e-03
   21    1518    2.636990e+03    5.303e-05    2.351e-04    5.481e-04
   22    1587    2.637047e+03    1.027e-05    1.885e-04    1.465e-04
   23    1656    2.637059e+03    0.000e+00    1.583e-04    3.424e-05
   24    1725    2.637059e+03    0.000e+00    2.416e-04    1.055e-06
   25    1794    2.637059e+03    0.000e+00    2.334e-04    2.272e-10
   26    1863    2.637059e+03    0.000e+00    6.141e-05    7.462e-15

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 111.767542 seconds.

Published with MATLAB® R2018a