【优化红绿灯】基于遗传算法优化绿灯时间实现单个十字路口的交通拥堵情况疏通附matlab代码

%% Starting point, clear everything in matlabtic;clear all;close all;clc;%% Problem FormulationFitnessFunction=@(C,g,x,c) TDi(C,g,x,c);     % FitnessFunctionnLights=4;                                   % Number of Traffic LightsnIntersections=1;                            % Number of Intersections (static as 1 intersection)VarSize=[1 nIntersections*nLights];          % Decision Chromosome genes based on number of IntersectionsgreenMin= 10;                                % Lower bound of GREEN LIGHTgreenMax= 60;                                % Upper bound of GREEN LIGHTCyclemin=60;                                 % Lower bound of CYCLECyclemax=180 ;RoadcapacityNSWE=[20,20,20,20];              % Road Capacity for NSWE respectivellyCarsNSWE=[20,20,11,17];RoadCongestion1NSWE=RoadcapacityNSWE-CarsNSWE;              % congestion according to free road spacesRoadCongestionNSWE=RoadCongestion1NSWE./RoadcapacityNSWE;   %  Volume/Capacity RATIOcarpass=5;%% Genetic Algorithm ParametersMaxIt=25;                                  % Maximum Number of IterationsnPop=400;                                     % Population Sizepc=0.5;                                      % Crossover Percentagenc=2*round(pc*nPop/2);                       % Number of Offsprings (parents)pm=0.02;                                      % Mutation Percentagenm=round(pm*nPop);                           % Number of Mutantsmu=0.1;                                      % Mutation Rate pinv=0.2;ninv=round(pinv*nPop);beta=8;                                      % Selection Pressure   %% Initialization% Individual Structureempty_individual.GreenNSWE=[];empty_individual.TotalDelay=[];% Population Structurepop=repmat(empty_individual,nPop,1);% Initialize Populationi=1;current_cycle=160-12; %estw kiklos 160 seconds - 12 seconds gia kitrinowhile i<=nPop          % Initialize Individual    pop(i).GreenNSWE=randi([greenMin greenMax],VarSize);        % Cycle time rules%     if(sum(CarsNSWE)<10)%         current_cycle(i)=randi([Cyclemin 80]);%     elseif(sum(CarsNSWE)<15)%         current_cycle(i)=randi([80 100]);%     elseif(sum(CarsNSWE)<20)%         current_cycle(i)=randi([100 120]);%     elseif(sum(CarsNSWE)<25)%         current_cycle(i)=randi([120 140]);%     elseif(sum(CarsNSWE)<30)%         current_cycle(i)=randi([140 160]);%     else%         current_cycle=180;%      end%     current_cycle=current_cycle(:);    if(sum(pop(i).GreenNSWE)>current_cycle)          continue;    end        % Individual Evaluation from Fitness Function    for j=1:nLights        % Measure Delay for each traffic light with current congestion        pop(i).TotalDelay(j)=FitnessFunction(current_cycle,pop(i).GreenNSWE(j),RoadCongestionNSWE(j),RoadcapacityNSWE(j));    end       % Summation of Total Delays quotients    pop(i).TotalDelay= real(sum(pop(i).TotalDelay));    i=i+1;end% Sort PopulationTotalDelay=[pop.TotalDelay];[TotalDelay, SortOrder]=sort(TotalDelay);pop=pop(SortOrder);% Store Best SolutionBestSol=pop(1);% Store Best FitnessBestDelay=pop(1).TotalDelay;% Worst FitnessWorstDelay=pop(end).TotalDelay;disp(['FIRST Population..........Best TotalDelay = ' num2str(BestDelay)]);    fprintf('\n')    disp('Green Timings in seconds:');    disp(['  North Green time = ' num2str(BestSol.GreenNSWE(1))]);    fprintf('\n')    disp(['  South Green time = ' num2str(BestSol.GreenNSWE(2))]);    fprintf('\n')    disp(['  West Green time = ' num2str(BestSol.GreenNSWE(3))]);    fprintf('\n')    disp(['  East Green time = ' num2str(BestSol.GreenNSWE(4))]);    fprintf('\n')%% Loop For Number of Iterationscount=0;for it=1:MaxIt       % TERMINATION CRITERIA IF NEEDED%     if(it>2)%         if(BestDelay(it-1)==BestDelay(it-2))%             count=count+1;%         else%         count=0;%         end%     end%     if(count>5)%         disp('5 Generations without evolution. Process Stopped !');%         break;%     end    % Calculate Selection Probabilities    P=exp(-beta*TotalDelay/WorstDelay);    P=P/sum(P);    %% Crossover    popc=repmat(empty_individual,nc/2,2);    k=1;    while k<=nc/2                % Select Parents Indices from roulette wheel        i1=RouletteWheelSelection(P);        i2=RouletteWheelSelection(P);                      % Select Parents        p1=pop(i1);        p2=pop(i2);         popc(k,1).GreenNSWE=p1.GreenNSWE;        popc(k,2).GreenNSWE=p2.GreenNSWE;        popc(k,1).TotalDelay=p1.TotalDelay;        popc(k,2).TotalDelay=p2.TotalDelay;        % Select random crossover point        i=randi([1 3]);                % crossover randomness         if(i==1)                              popc1=popc(k,1).GreenNSWE(2);                popc(k,1).GreenNSWE(2)= popc(k,2).GreenNSWE(2);                popc(k,2).GreenNSWE(2)=popc1;                            popc1=popc(k,1).GreenNSWE(3);                popc(k,1).GreenNSWE(3)= popc(k,2).GreenNSWE(3);                popc(k,2).GreenNSWE(3)=popc1;                        popc1=popc(k,1).GreenNSWE(4);                popc(k,1).GreenNSWE(4)= popc(k,2).GreenNSWE(4);                popc(k,2).GreenNSWE(4)=popc1;          elseif(i==2)                    popc1=popc(k,1).GreenNSWE(3);                popc(k,1).GreenNSWE(3)= popc(k,2).GreenNSWE(3);                popc(k,2).GreenNSWE(3)=popc1;                        popc1=popc(k,1).GreenNSWE(4);                popc(k,1).GreenNSWE(4)= popc(k,2).GreenNSWE(4);                popc(k,2).GreenNSWE(4)=popc1;           else                             popc1=popc(k,1).GreenNSWE(4);                popc(k,1).GreenNSWE(4)= popc(k,2).GreenNSWE(4);                popc(k,2).GreenNSWE(4)=popc1;        end                % check if new green times are out constraints 10-60s. If it is        % get it to closer min or max        popc(k,1).GreenNSWE=max(popc(k,1).GreenNSWE,greenMin);        popc(k,1).GreenNSWE=min(popc(k,1).GreenNSWE,greenMax);        popc(k,2).GreenNSWE=max(popc(k,2).GreenNSWE,greenMin);        popc(k,2).GreenNSWE=min(popc(k,2).GreenNSWE,greenMax);                if(sum(popc(k,1).GreenNSWE)>current_cycle || sum(popc(k,2).GreenNSWE)>current_cycle)            continue;        end                % Evaluate Generated Offsprings for each traffic light according to        % the corresponding traffic congestion        for j=1:nLights            popc(k,1).TotalDelay(j)=FitnessFunction(current_cycle, popc(k,1).GreenNSWE(j), RoadCongestionNSWE(j),RoadcapacityNSWE(j));            popc(k,2).TotalDelay(j)=FitnessFunction(current_cycle, popc(k,2).GreenNSWE(j), RoadCongestionNSWE(j),RoadcapacityNSWE(j));        end                % TOTAL DELAY which correspongs to the summation of of 4 lights quotients        popc(k,1).TotalDelay= real(sum(popc(k,1).TotalDelay));        popc(k,2).TotalDelay= real(sum(popc(k,2).TotalDelay));             k=k+1; %step    end        % Make 2 rows 1    popc=popc(:);    % Sort popc matrix according to TotalDelay    TotalDelay=[popc.TotalDelay];    [TotalDelay, SortOrder]=sort(TotalDelay);    popc=popc(SortOrder);        %% Mutation    % Create empty Matrix with length the number of mutants     popm=repmat(empty_individual,nm,1);    k=1;    while k<=nm        % Select Parent population        i=randi([1 nPop]); %nPop value 100        p=pop(i);                % Apply Mutation           nVar=4;        nmu=ceil(mu*nVar);        j=randi([1 nVar]);        prosimo=randi([-1 1]);        sigma=prosimo*0.02*(greenMax-greenMin);                mutated=p.GreenNSWE(j)+sigma;                popm(k).GreenNSWE = p.GreenNSWE;        popm(k).GreenNSWE(j)=mutated;                popm(k).GreenNSWE=max(popm(k).GreenNSWE,greenMin);        popm(k).GreenNSWE=min(popm(k).GreenNSWE,greenMax);                if(sum(popm(k).GreenNSWE)>current_cycle)            continue;        end                 for j=1:nLights            % Evaluate Mutant             popm(k).TotalDelay(j)=FitnessFunction(current_cycle, popm(k).GreenNSWE(j), RoadCongestionNSWE(j),RoadcapacityNSWE(j));        end        % Summation of delay quotients        popm(k).TotalDelay=real(sum(popm(k).TotalDelay));               k=k+1; %step    end            %% INVERSION    % Create empty Matrix     popinv=repmat(empty_individual,nm,1);    k=1;    while k<=ninv                % Select Parent population        i=randi([1 nPop]);        p=pop(i);                % Apply Inversion                nVar=numel(p.GreenNSWE);        randomgene1=randi([1 4]);        randomgene2=randi([1 4]);        y=p.GreenNSWE;        popinv(k).GreenNSWE=y;        x=popinv(k).GreenNSWE(randomgene1);        popinv(k).GreenNSWE(randomgene1)=popinv(k).GreenNSWE(randomgene2);              popinv(k).GreenNSWE(randomgene2)=x;        popinv(k).GreenNSWE=max(popinv(k).GreenNSWE,greenMin);        popinv(k).GreenNSWE=min(popinv(k).GreenNSWE,greenMax);        if(sum(popinv(k).GreenNSWE)>current_cycle)            continue;        end                for j=1:nLights            % Evaluate Mutant             popinv(k).TotalDelay(j)=FitnessFunction(current_cycle, popinv(k).GreenNSWE(j), RoadCongestionNSWE(j),RoadcapacityNSWE(j));        end        % Summation of delay quotients        popinv(k).TotalDelay=real(sum(popinv(k).TotalDelay));                      k=k+1;    end    % Make 2 rows 1    popinv=popinv(:);    %% Merge Population    pop=[pop        popc        popm        popinv]; %#ok         % Sort New Population according to TotalDelay    TotalDelay=[pop.TotalDelay];    [TotalDelay, SortOrder]=sort(TotalDelay);    pop=pop(SortOrder);        % Update Worst Cost    WorstDelay=max(WorstDelay,pop(end).TotalDelay);        % Keep the Best Population from the given number    pop=pop(1:nPop);    TotalDelay=TotalDelay(1:nPop);        % Store Best Solution Ever Found    BestSol=pop(1);        % Store Best Cost Ever Found    BestDelay(it)=BestSol.TotalDelay;            % Show Iteration Information    disp(['                   Iteration ' num2str(it) ': Best TotalDelay = ' num2str(BestDelay(it))]);    fprintf('\n')    disp('Green Timings:');    fprintf('\n')    disp(['  North Green time = ' num2str(BestSol.GreenNSWE(1))'' ' seconds']);    fprintf('\n')    disp(['  South Green time = ' num2str(BestSol.GreenNSWE(2))'' ' seconds']);    fprintf('\n')    disp(['  West Green time = ' num2str(BestSol.GreenNSWE(3))'' ' seconds']);    fprintf('\n')    disp(['  East Green time = ' num2str(BestSol.GreenNSWE(4))'' ' seconds']);    fprintf('\n')    %end of generation    end    disp(' ****************************************************************' );    disp('  CASE: Every 5 seconds 2 vehicles leaves the corresponding road ' );    disp('  Expected vehicles left through North road' );    disp(round(2*BestSol.GreenNSWE(1)/carpass));    disp('  Expected vehicles left through South road' );    disp(round(2*BestSol.GreenNSWE(2)/carpass));    disp('  Expected vehicles left through West road' );    disp(round(2*BestSol.GreenNSWE(3)/carpass));    disp('  Expected vehicles left through East road' );    disp(round(2*BestSol.GreenNSWE(4)/carpass));    fprintf('\n')    disp(' ****************************************************************' );disp(['Cycle Time = ' num2str(current_cycle)'' ' seconds']);%% Results / Plotsfigure(1);semilogy(BestDelay,'LineWidth',2);% plot(BestCost,'LineWidth',2);xlabel('Iteration');ylabel('Total Delay');grid on;toc