GAconstrain(prob)

function r = GAconstrain(prob)
% Genetic Algorithm(real coding)
% By: Javad Ivakpour
% E-mail: [e-mail]
% May 2006

r = struct('alg', 'Genetic Algorithm(real coding)',...
    'authors', 'Javad Ivakpour [e-mail]',...
    'lastChanger', 'Dmitrey, [e-mail]');

%------------------------        parameters        ------------------------
var=prob.n;            % Number of variables (this item must be equal to the
                  %   number of variables that is used in the function in
                  %   fun00.m file)
n=100;            % Number of population

m0=30;            %  Number of generations that max value remains constant
                  %   (use for termination criteria)
nmutationG=20;                  %number of mutation children(Gaussian)
nmutationR=20;                  %number of mutation children(random)
nelit=2;                        %number of elitism children
valuemin=prob.lb;     % min possible value of variables
valuemax=prob.ub;      % max possible value of variables
if any(isinf(prob.lb)) || any(isinf(prob.ub))
    prob.err('solver GAconstrain requires prob.lb & prob.ub')
end
A=prob.A;     % constrains: Ax+B must be less or equal to 0
B=-prob.b;
%-------------------------------------------------------------------------

%checking are any solver parameters is structure prob.(name of solver) 
ExtractRoutineParamsFromProb;%script-file
%if present - they will replace defaults 

nmutation=nmutationG+nmutationR;
sigma=(valuemax-valuemin)/10;    %Parameter that related to Gaussian
                                 %   function and used in mutation step
max1=zeros(nelit,var);
parent=zeros(n,var);

p = zeros(n,var);
for l=1:var
    p(:,l)=valuemin(l)+rand(n,1).*(valuemax(l)-valuemin(l));
end
if isreal(prob.x0); p(1,:) = prob.x0'; end

m=m0;
maxvalue=ones(m,1)*-1e10;
maxvalue(m)=-1e5;
g=0;

niter = 0;
no_feasible = -realmax/1e2;
%-------------   ****    termination criteria   ****-------------
% while abs(maxvalue(m)-maxvalue(m-(m0-1)))>0.001*maxvalue(m) &&...
%         (abs(maxvalue(m))>1e-10 && abs(maxvalue(m-(m0-1)))>1e-10)...
%         && m<10000 && abs(maxvalue(m)-meanvalue(m))>1e-5 || niter<prob.maxIter
while 1 % by Dmitrey
    sigma=sigma./(1.05);% reducing the sigma value
    %  ------  **** % reducing the number of mutation()random   **** ----
    g=g+1;
    if g>10 && nmutationR>0
        g=0;
        nmutationR=nmutationR-1;
        nmutation=nmutationG+nmutationR;
    end


    %-------------   ****    function evaluation   ****-------------
    for i=1:n%CHANGED BY DMITREY
        x = p(i,:);
        x = x(:);
        nfz = 0;
        if length(prob.b)
            for j=1:length(prob.b)
                if any(A(j,:)*x+B(j)>0)
                    y(i)=no_feasible;
                    nfz = 1;
                    break
                end
            end
        end
        if ~nfz && isfield(prob, 'c') && ~isempty(prob.c) && any(prob.c(x, prob)>prob.TolCon)
            y(i)=no_feasible;
            nfz = 1; 
        end
        if ~nfz; y(i) = -prob.f(x(:), prob); end
    end
    ind_n = y==no_feasible;
    if all(ind_n==1)
        r.istop = -11;
        return
    end
    miny = min(y);
    y(ind_n) = miny - 1e-8 * abs(miny);

    s=sort(y);
    maxvalue1(1:nelit)=s(n:-1:n-nelit+1);
    if nelit==0
        maxvalue1(1)=s(n);
        for i=1:n
            if y(i)==maxvalue1(1)
                max1(1,:)=p(i,:);
            end
        end
    end
    for k=1:nelit
        for i=1:n
            if y(i)==maxvalue1(k)
                max1(k,:)=p(i,:);
            end
        end
    end
    y=y-miny;
    sy = sum(y);
    sumd=y./sy;



    %-------------   ****   Selection: Rolette wheel   ****-------------
    for l=1:n
        sel=rand;
        sumds=0;
        j=1;
        while sumds<sel
            sumds=sumds+sumd(j);
            j=j+1;
        end
        parent(l,:)=p(j-1,:);
    end
    p=zeros(n,var);

    %-------------   ****    regeneration   ****-------------
    for l=1:var


        %-------------   ****    cross-over   ****-------------
        for j=1:ceil((n-nmutation-nelit)/2)
            t=rand*1.5-0.25;
            p(2*j-1,l)=t*parent(2*j-1,l)+(1-t)*parent(2*j,l);
            p(2*j,l)=t*parent(2*j,l)+(1-t)*parent(2*j-1,l);
        end


        %-------------   ****    elitism   ****-------------
        for k=1:nelit
            p((n-nmutation-k+1),l)=max1(k,l);
        end

        %------------****mutation by Dmitrey****-------------
        phi=1-2*rand(nmutation-nmutationR,1);
        z = erfinv(phi)*sqrt(2);
        p(n-nmutation+1:n-nmutationR,l) = z * sigma(l)+parent(n-nmutation+1:n-nmutationR,l);
% 
%         %-------------   ****    mutation   ****-------------
%         for i=n-nmutation+1:n-nmutationR
%             phi=1-2*rand;
%             z=erfinv(phi)*(2^0.5);
%             p(i,l)=z*sigma(l)+parent(i,l);
% 
%         end
        for i=n-nmutationR+1:n
            p(i,1:var)=(valuemin(1:var)+rand(var,1).*(valuemax(1:var)...
                -valuemin(1:var)))';
        end
        for i=1:n
            %for l=1:var %commented by Dmitrey
            %cycle by l already exists
            if p(i,l)<valuemin(l)
                p(i,l)=valuemin(l);
            elseif p(i,l)>valuemax(l)
                p(i,l)=valuemax(l);
            end
            %end  %commented by Dmitrey
        end
    end
    m=m+1;
    niter = niter+1;
    maxvalue(m)=maxvalue1(1);


    r.f = -maxvalue1(1);
    r.x = max1(1,:);
    r.x = r.x(:);
    r = prob.iterfcn(prob, r);
    if r.istop; return; end
end

% clc
% disp('     Genetic Algorithm(real coding)   ')
% disp('          By: Javad Ivakpour          ')
% disp('       E-mail: [e-mail]       ')
% disp('**************************************')
% num_of_fun_evaluation=n*m