testfunction11.m

Colonial Competitive Algorithm

% Test functions for optimization
% These are the test functions that appear in Appendix I.
% Set funnum to the function you want to use.
% funnum=17 is for a MOO function
% Haupt & Haupt
% 2003
function f=testfunction11(x)  % 20

% 22 MIMO PID Tunuing

bias = 2;

funnum = -1;
alpha = 2;
if funnum==-1 %F1 Rasrigin
    f = 18.3549 + 20 + x(:,1).^2 + x(:,2).^2 - 10 * ( cos(alpha*pi*x(:,1)) + cos(alpha*pi*x(:,2)) );
elseif funnum==0 %F1
    f = (x(:,1) - bias).^2 + (x(:,2) - bias).^2;
elseif funnum==1 %F2
    f=abs(x)+cos(x);

elseif funnum==2 %F2
    f=abs(x)+sin(x);
elseif funnum==3 %F3
    f=x(:,1).^2+x(:,2).^2;
elseif funnum==4 %F4
    f=100*(x(:,2).^2-x(:,1)).^2+(1-x(:,1)).^2;
elseif funnum==5 %F5
    f(:,1)=sum(abs(x')-10*cos(sqrt(abs(10*x'))))';
elseif funnum==6 %F6
    f=(x.^2+x).*cos(x);
elseif funnum==7 %F7
    %     x = x * 10;
    f=x(:,1).*sin(4*x(:,1))+1.1*x(:,2).*sin(2*x(:,2));
elseif funnum==8 %F8
    f=x(:,2).*sin(4*x(:,1))+1.1*x(:,1).*sin(2*x(:,2));
elseif funnum==9 %F9
    f(:,1)=x(:,1).^4+2*x(:,2).^4+randn(length(x(:,1)),1);
elseif funnum==10 %F10
    f(:,1)=20+sum(x'.^2-10*cos(2*pi*x'))';
elseif funnum==11 %F11
    f(:,1)=1+sum(abs(x').^2/4000)'-prod(cos(x'))';
elseif funnum==12 %F12
    f(:,1)=.5+(sin(sqrt(x(:,1).^2+x(:,2).^2).^2)-.5)./(1+.1*(x(:,1).^2+x(:,2).^2));
elseif funnum==13 %F13
    aa=x(:,1).^2+x(:,2).^2;
    bb=((x(:,1)+.5).^2+x(:,2).^2).^0.1;
    f(:,1)=aa.^0.25.*sin(30*bb).^2+abs(x(:,1))+abs(x(:,2));
elseif funnum==14 %F14
    f(:,1)=besselj(0,x(:,1).^2+x(:,2).^2)+abs(1-x(:,1))/10+abs(1-x(:,2))/10;
elseif funnum==15 %F15
    f(:,1)=-exp(.2*sqrt((x(:,1)-1).^2+(x(:,2)-1).^2)+(cos(2*x(:,1))+sin(2*x(:,1))));
elseif funnum==16 %F16
    f(:,1)=x(:,1).*sin(sqrt(abs(x(:,1)-(x(:,2)+9))))-(x(:,2)+9).*sin(sqrt(abs(x(:,2)+0.5*x(:,1)+9)));
elseif funnum==17 %MOO function
    x=x+1;
    f(:,1)=(x(:,1)+x(:,2).^2+sqrt(x(:,3))+1./x(:,4))/8.5;
    f(:,2)=(1./x(:,1)+1./x(:,2)+x(:,3)+x(:,4))/6;
elseif funnum == 18
%     xN = (x-0.5)*2;
%     xN = (x);
x = (x)*12;
    x1 = x(:,1);
    x2 = x(:,2);
    
    x1Min = 3;
    x2Min = 3;
    f = 20 + (x1 - x1Min).^2 + (x2 - x2Min).^2 - 10 * (cos(2*pi*(x1 - x1Min)) + cos(2*pi*(x2 - x2Min)));
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif funnum == 20 % Game Theory
    A = [20     20      20
        25     25      25
        30     30      30];
    for i = 1 : size(x,1)
        f(i,1) = x(i,:)*A(:,3) - x(i,:)*A*x(i,:)';
    end
    f=x(:,1).^2+x(:,2).^2;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
elseif funnum == 21 % PID Tuning
    s = tf('s');
    %     G = tf([1 2],[1 5 5 5]);
    G = 4.228/(s+0.5)/(s^2+1.64*s+8.456);
    for jj = 1:size(x,1)

        KP = x(jj,1);
        KI = x(jj,2);
        KD = x(jj,3);
        PID_Controller = KP + KD*s + KI/s;
        % G = (1 - 5*s)  /  (  (1+10*s)*(1+20*s)  )
        %
        % Arash System
        SERIES = series(G,PID_Controller);
        ClosedLoop = feedback(SERIES,1);
        gam = 0.01;
        t_in = [0:gam:40];
        [y,t] = step(ClosedLoop,t_in);
        %         step(ClosedLoop,t_in);
        % hold on
        % step(G);
        % [y,t] = step(G,t_in);

        %=====================================================
        %  Finding Max Overshoot
        %=====================================================
        maxovershoot = 0;
        if max(y)>1
            ind = find(y==max(y));
            maxovershoot = (y(ind) - 1);
            %     hold on
            %     plot(t(ind),y(ind),'rp')
            %     hold off
        end

        %=====================================================
        %  Finding Rise Time
        %=====================================================
        y_1_10 = abs(y - 0.1);
        y_9_10 = abs(y - 0.9);

        index1_10 = find(   y_1_10 < (5*min(y_1_10))  );
        index1_10 = index1_10(1);

        index9_10 = find(   y_9_10 < (5*min(y_9_10))  );
        index9_10 = index9_10(1);

        RiseTime = t(index9_10)-t(index1_10);y_1_10 = abs(y - 0.1);
        % hold on
        % plot(t(index9_10),y(index9_10),'gp' , t(index1_10),y(index1_10),'gp')
        % hold off

        %=====================================================
        %  Finding Settling Time
        %=====================================================
        y_1_2_100 = abs(y - 1.05);
        y_0_98_100 = abs(y - 0.95);
        y_1 = abs(y - 1);

        index1_2_100 = find(   y_1_2_100 < (0.01 + min(y_1_2_100))  );
        index0_98_100 = find(   y_0_98_100 < (0.01 + min(y_0_98_100))  );
        % hold on
        % plot(t(index1_2_100),y(index1_2_100),'k.'  ,  t(index0_98_100),y(index0_98_100),'y.')
        % hold off
        settling_time = t(end);
        settling_time_index = numel(t);
        all_indexes = sort([index1_2_100 ;  index0_98_100]);

        for ii = 1:numel(all_indexes)
            if (          sum( (y_1(all_indexes(ii):end)>0.05) )  ) == 0
                settling_time_index = all_indexes(ii);
                settling_time = t(all_indexes(ii));
                break
            end
        end
        %         hold on
        %         plot(settling_time,y(settling_time_index),'kp')
        %         hold off

        %=====================================================
        %  Finding IAE
        %=====================================================
        IAE = sum(abs(y-1))*gam;
        %=====================================================
        %  Returning
        %=====================================================
        f(jj,1) = 1*(maxovershoot  +  RiseTime +  IAE  +  settling_time);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MIMO PID Tuning %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

elseif funnum == 22 % MIMO PI Tuning
    %     x = (x-.5) * 2;
    s = tf('s');
    %   PADE APPROXIMATION IS USED TO APROXIMATE EXP(-TS)
    %         G11 = 12.8*((1-.5*s+.5*(-.5*s)^2)/(1+.5*s+.5*(.5*s)^2))/(16.7*s+1);
    %         G12 = -18.9*((1-1.5*s+.5*(-1.5*s)^2)/(1+1.5*s+.5*(1.5*s)^2))/(21*s+1);
    %         G21 = 6.6*((1-3.5*s+.5*(-3.5*s)^2)/(1+3.5*s+.5*(+3.5*s)^2))/(10.9*s+1);
    %         G22 = -19.4*((1-1.5*s+.5*(-1.5*s)^2)/(1+1.5*s+.5*(1.5*s)^2))/(14.4*s+1);
    %
    ApproxOrder = 2;
    [n1,d1] = pade(1,ApproxOrder);
    [n3,d3] = pade(3,ApproxOrder);
    [n7,d7] = pade(7,ApproxOrder);

    G11 = 12.8*tf(n1,d1)/(16.7*s+1);
    G12 = -18.9*tf(n3,d3)/(21*s+1);
    G21 = 6.6*tf(n7,d7)/(10.9*s+1);
    G22 = -19.4*tf(n3,d3)/(14.4*s+1);

    G = [ G11 G12;  G21 G22];
    gam  =  .1;
    t_in = [0:gam:800];

    KP11 = x(:,1);
    KI11 = x(:,2);
    KP12 = x(:,3);
    KI12 = x(:,4);
    KD12 = x(:,5);
    KP21 = x(:,6);
    KI21 = x(:,7);
    KD21 = x(:,8);
    KP22 = x(:,9);
    KI22 = x(:,10);


    for jj = 1:size(x,1)
        jj
        %         x = [KP11 KI11 | KP12 KI12 KD12 | KP21 KI21 KD21 | KP22 KI22 ]
        %
        %         KP11 = 0.184; KI11 = 0.184/3.92; KP12 = -0.0102; KI12 = -0.0102/0.445; KD12 = -0.0102*(-0.804);...
        %         KP21 = -0.0674; KI21 = -0.0674/(-4.23); KD21 = -0.0674 * 0.796; KP22 = -0.0660; KI22 = -0.0660/4.25; jj = 1; %Auto Tuning Controller
        % % %
        % %         KP11 = 0.9971;        KI11 = 0.0031;        KP22 = -.0141;        KI22 = -0.0071; % Their proposed algorithm
        %
        %        KP11 = 0.5511;        KI11 = 0.0018;        KP22 = -.0182;        KI22 = -0.0067; % Traditional GA

        %         KP11 = 0.3750;        KI11 = 0.0452;        KP22 = -.0750;        KI22 = -0.0032; % BLT Method
        Controller11 = KP11(jj) + KI11(jj)/s;
        Controller12 = KP12(jj) + KI12(jj)/s + KD12(jj)*s;
        Controller21 = KP21(jj) + KI21(jj)/s + KD21(jj)*s;
        Controller22 = KP22(jj) + KI22(jj)/s;

        K = [ Controller11 Controller12;  Controller21 Controller22];

        TF = G*K*inv(eye(2,2)+G*K);

        [Y,t] = step(TF,t_in);
        y11 = Y(:,1,1);
        y12 = Y(:,1,2);
        y21 = Y(:,2,1);
        y22 = Y(:,2,2);
        %         step(TF,t_in);
        %=====================================================
        %  Finding IAE
        %=====================================================
        IAE = sum(abs(y11-1) + abs(y12-0) + abs(y21-0) + abs(y22-1)) * gam
        %=====================================================