gadelphi.pas

关于加气站的遗传算法实现

unit GADelphi;
interface
uses
  Math, SysUtils,Classes;

type

  PUnsigned = array of Cardinal;

  Pint = ^integer;
  individual = record
    chrom: PUnsigned;
    fitness: double;
    varible: double;
    xsite: integer;
    parent: array[0..1] of integer;
    utility: Pint;
  end;
 //////2005-4-1 Add---------------------
  Genes = record  //一个基因，也就是表示一个需求点
        x:single;
        y:single;
        Demand: single; //这个点的需求量
        DPorSP :integer;   //改点为需求点或者供应点，0为需求点，1为供应点
        Supplypoint: string; //该点如果为供应点，它所供应的基因的下标字符串
        Distance: single;   //该点到供应点的距离，若该点为供应点则为0
    end;
  //////2005-4-1 Add---------------------
  bestever = record
    chrom: PUnsigned;
    fitness: Double;
    varible: Double;
    generation: Integer;
  end;
  //* 随机数发生器使用的静态变量 */

  TGA = class(Tobject)
  public
    oldpop, temp, newpop: array of individual;
    //////2005-4-1 Add---------------------
    g,d,p,V,k,Q : single;//汽车的时间和距离费用参数
    MyGenes : array of Genes;
    FData,FResult: TextFile ;//数据输入输出文件
    supplypoint: array of integer;
    //////2005-4-1 Add---------------------
    bestfit: bestever;
    sumfitness: double;
    max, avg, min: double;
    pcross, pmutation: single;
    popsize, lchrom, chromsize, gen, maxgen, run, maxruns: integer;
    printstrings: integer; //* 输出染色体编码的判断，0 -- 不输出, 1 -- 输出 */
    nmutation: integer; //* 当前代变异发生次数 */
    ncross: integer;

    procedure initpop;
    procedure objfunc(var critter: individual);
    function flip(prob: single): boolean;
    function select: integer;
    function crossover(var parent1, parent2, child1, child2: PUnsigned): integer;
    function rnd(low, high: integer): integer;
    procedure mutation(child: PUnsigned);
    procedure generation;
    procedure preselect;
    constructor Create;
    procedure initdata;
    procedure initmalloc;
    procedure statistics(pop: array of individual);
    procedure initreport;
    procedure skip(skipcount: integer);
    procedure Go;
    procedure report;
    procedure repchar(ch: pchar; repcount: integer);
    procedure writepop;
    procedure writechrom(chrom: PUnsigned);
    procedure GARandomize;
    procedure advance_random;
    function randomperc: single;
    function randomnormaldeviate: double;
    procedure warmup_random(random_seed: single);
    function SplitString(const Source,ch:string):integer;////
    procedure reportsupplypoint(var critter: PUnsigned);
  end;

implementation
var
  oldrand: array[0..54] of double;
  jrand: integer;
  rndx2: double;
  rndcalcflag: integer;

procedure TGA.GARandomize; // * 设定随机数种子并初始化随机数发生器 * /
var
  randomseed: single;
  j1: integer;
begin
  for j1 := 0 to 54 do
    oldrand[j1] := 0.0;
  jrand := 0;
  //repeat
    //write('随机数种子[0 - 1]: ');
    //readln(randomseed);
  //until not ((randomseed < 0.0) and (randomseed > 1.0));
  randomseed := 0.5;
  warmup_random(randomseed);
end;

procedure TGA.warmup_random(random_seed: single); // * 初始化随机数发生器 * /
var
  j1, ii: integer;
  new_random, prev_random: double;
begin
  oldrand[54] := random_seed;
  new_random := 0.000000001;
  prev_random := random_seed;
  for j1 := 1 to 54 do
  begin
    ii := (21 * j1) mod 54;
    oldrand[ii] := new_random;
    new_random := prev_random - new_random;
    if (new_random < 0.0) then
      new_random := new_random + 1.0;
    prev_random := oldrand[ii];
  end;
  advance_random;
  advance_random();
  advance_random();
  jrand := 0;
end;

procedure TGA.advance_random; //* 产生55个随机数 */
var
  j1: integer;
  new_random: double;
begin
  for j1 := 0 to 23 do
  begin
    new_random := oldrand[j1] - oldrand[j1 + 31];
    if (new_random < 0.0) then
      new_random := new_random + 1.0;
    oldrand[j1] := new_random;
  end;
  for j1 := 24 to 54 do
  begin
    new_random := oldrand[j1] - oldrand[j1 - 24];
    if (new_random < 0.0) then new_random := new_random + 1.0;
    oldrand[j1] := new_random;
  end;
end;

function TGA.randomnormaldeviate: double; //* 产生随机标准差 */
var
  t, rndx1: double;
begin

  if rndcalcflag = 1 then
  begin
    rndx1 := sqrt(-2.0 * log2(randomperc));
    t := 6.2831853072 * randomperc;
    rndx2 := rndx1 * sin(t);
    rndcalcflag := 0;
    result := (rndx1 * cos(t));
  end
  else
  begin
    rndcalcflag := 1;
    result := (rndx2);
  end;
end;

function TGA.randomperc: single;
  //*与库函数random()作用相同, 产生[0,1]之间一个随机数 */
begin
  inc(jrand, 1);
  if (jrand >= 55) then
  begin
    jrand := 1;
    advance_random();
  end;
  result := oldrand[jrand];
end;

constructor TGA.Create; //* 遗传算法初始化 */
begin
  //初始化输出文件；
   AssignFile(FResult, 'result.txt');
   rewrite(FResult);

     //* 键盘输入遗传算法参数 */
  initdata;
    //* 确定染色体的字节长度 */
  chromsize := (lchrom div (8 * sizeof(cardinal)));
  if ((lchrom mod (8 * sizeof(cardinal))) = 0) or (lchrom < (8 *
    sizeof(cardinal)))  then
    inc(chromsize, 1);
    //*分配给全局数据结构空间 */
  initmalloc;
    //* 初始化随机数发生器 */
  GArandomize;
    //* 初始化全局计数变量和一些数值*/
  nmutation := 0;
  ncross := 0;
  bestfit.fitness := -1E10;
  bestfit.generation := 0;
    //* 初始化种群，并统计计算结果 */
  initpop;
  statistics(oldpop);
  report;
  initreport;
end;

procedure TGA.initreport; //* 初始参数输出 /*
begin
  skip(1);
  writeln(FResult, '             基本遗传算法参数\n');
  writeln(FResult, ' -------------------------------------------------');
  writeln(FResult, format(' 种群大小(popsize) = %d', [popsize]));
  writeln(FResult, format(' 染色体长度(lchrom) = %d', [lchrom]));
  writeln(FResult, format(' 最大进化代数(maxgen) = %d', [maxgen]));
  writeln(FResult, format(' 交叉概率(pcross) = %8.6f', [pcross]));
  writeln(FResult, format(' 变异概率(pmutation) = %8.6f', [pmutation]));
  writeln(FResult, ' -------------------------------------------------');
  skip(1);

end;

procedure TGA.skip(skipcount: integer);
var j: integer;
begin
  for j := 1 to skipcount do
    writeln(FResult,'');
end;

procedure TGA.statistics(pop: array of individual); //* 计算种群统计数据 */
var
  i, j: integer;
begin
  sumfitness := 0.0;
  min := pop[0].fitness;
  max := pop[0].fitness;
    //* 计算最大、最小和累计适应度 */
  for j := 0 to popsize - 1 do
  begin
    sumfitness := sumfitness + pop[j].fitness;
    if (pop[j].fitness > max) then
      max := pop[j].fitness;
    if (pop[j].fitness < min) then
      min := pop[j].fitness;
    // * new global best - fit individual * /
    if (pop[j].fitness > bestfit.fitness) then
    begin
      for i := 0 to chromsize - 1 do
        bestfit.chrom[i] := pop[j].chrom[i];
      bestfit.fitness := pop[j].fitness;
      bestfit.varible := pop[j].varible;
      bestfit.generation := gen;
    end;
  end;
  // * 计算平均适应度 * /
  avg := sumfitness / popsize;
end;

procedure TGA.initmalloc;
var
  nbytes: cardinal;
  j: integer;
begin
//  char  *malloc();
  //* 分配给当前代和新一代种群内存空间 */
 // nbytes = popsize * sizeof(individual);
{  if((oldpop = (struct individual *) malloc(nbytes)) == NULL)
    nomemory("oldpop")}
  setlength(oldpop, popsize);
  setlength(newpop, popsize);

{  if((newpop = (struct individual *) malloc(nbytes)) == NULL)
    nomemory("newpop");}
  //* 分配给染色体内存空间 */
//  nbytes = chromsize * sizeof(unsigned);
  for j := 0 to popsize - 1 do
  begin
    {  if((oldpop[j].chrom = (unsigned *) malloc(nbytes)) == NULL)
 nomemory("oldpop chromosomes");}
    SetLength(oldpop[j].chrom, chromsize);
    SetLength(newpop[j].chrom, chromsize);
    {  if((newpop[j].chrom = (unsigned *) malloc(nbytes)) == NULL)
 nomemory("newpop chromosomes");}
  end;
  SetLength(bestfit.chrom, chromsize);
end;

procedure TGA.initdata; //* 遗传算法参数输入 */
var
  answer: array[0..1] of char;
  i:integer;
  tmp1,tmp2,tmp3:single;
begin
  //////2005-4-1 Add---------------------
   //chu shi hua shu ju wenjian
   AssignFile(FData, 'data.txt');
   reset(FData);
  //////2005-4-1 Add---------------------
  //writeln('种群大小(20-100)：');
  readln(FData,popsize);
  if ((popsize mod 2) <> 0) then
  begin
    writeln('种群大小已设置为偶数');
    inc(popsize, 1);
  end;
  //writeln('染色体长度(8-40)：');
  readln(FData,lchrom);
  //writeln('是否输出染色体编码(y/n)');
  printstrings := 1;
  //readln(answer);
  //if (comparetext(answer, 'n') = 0) then
   // printstrings := 0;
  //writeln('最大世代数(100 - 300)：');
  readln(FData,maxgen);
  //writeln('交叉率(0.2 - 0.9)：');
  readln(FData,pcross);
  //writeln('变异率(0.01 - 0.1)：');
  readln(FData,pmutation);
//////2005-4-1 Add---------------------

//此处输入汽车时间和距离费用参数
// g表示CNG汽车百公里耗气量（已知）;
       readln(FData,g);
//p表示CNG的单位价格（已知）;
        readln(FData,p);
//V表示CNG汽车前往加气的平均行驶速度（已知）;
        readln(FData,v);
//k表示CNG汽车的时间价值（已知）;
        readln(FData,k);
//Q表示CNG汽车单车标准加气量（已知）；
        readln(FData,Q);
//此处增加坐标数据的输入，通过文件输入；
setlength(MyGenes, lChrom);
setlength(supplypoint, lchrom);

  for i:=0 to lchrom-1 do
  begin
      Readln(FData,tmp1, tmp2, tmp3);
      MyGenes[i].x := tmp1;
      MyGenes[i].y := tmp2;
      MyGenes[i].Demand := tmp3;
  end;
closefile(FData);

//////2005-4-1 Add---------------------
end;

procedure TGA.initpop;
var
  j, j1, k, stop: integer;
  mask: Cardinal;
begin
  mask := 1;
  for j := 0 to popsize - 1 do
  begin
    for k := 0 to chromsize - 1 do
    begin
      oldpop[j].chrom[k] := 0;
      if k = chromsize - 1 then
        stop := lchrom - (k * (8 * sizeof(cardinal)))
      else
        stop := 8 * sizeof(cardinal);
      for j1 := 1 to stop do
      begin
        oldpop[j].chrom[k] := oldpop[j].chrom[k] shl 1;
        if flip(0.5) then
          oldpop[j].chrom[k] := oldpop[j].chrom[k] or mask;
      end;
      oldpop[j].parent[0] := 0;
      oldpop[j].parent[1] := 0;
      oldpop[j].xsite := 0;
      objfunc(oldpop[j]);
    end;
  end;
end;

procedure TGA.objfunc(var critter: individual); //* 计算适应度函数值 */
var
  mask, bitpos, tp: cardinal;
  bitpow: double;
  i,j, k, stop: integer;

//////2005-4-1 Add---------------------
  DPArray: array of integer;  //一个染色体中需求点的下标数组
  SPArray: array of integer; //一个染色体中供应点的下标数组
  DPNum,SPNum:  integer;//  一个染色体中需求点和供应点 的个数
  tmpDistance,tmpx,tmpy,tmpDistance2: single;//
  sp:integer; //某一需求点的供应点的下标
  StationDemand: single;
  StationCost,totalStationCost: single;//建站成本；
  TimeAndDisCost: single;
  totalDemand,totalSupply: double;
begin
  DPNum :=0;
  SPNum :=0;
  totalDemand:=0;
  totalSupply:=0;
  setlength(DPArray,lchrom);
  setlength(SPArray,lchrom);

  /////2005-4-1 Add---------------------
  mask := 1;
  critter.varible := 0.0;
  for k:=0 to lchrom-1 do
  begin
      mygenes[k].DPorSP :=0;
      mygenes[k].Distance :=0;
      mygenes[k].Supplypoint := '';
      totalDemand := totalDemand+ mygenes[k].Demand;
  end;
  for k := 0 to chromsize - 1 do
  begin
    if (k = (chromsize - 1)) then
      stop := lchrom - (k * (8 * sizeof(cardinal)))
    else
      stop := 8 * sizeof(cardinal);
    tp := critter.chrom[k];
    for j := 0 to stop - 1 do
    begin
      bitpos := j + (8 * sizeof(cardinal)) * k;
      if ((tp and mask) = 1) then
      begin
   //////2005-4-1 Add---------------------
        Mygenes[bitpos].DPorSP := 1;
        SPArray[SPNum] := bitpos ;
        SPNum := SPnum +1
      end
      else
      begin
        Mygenes[bitpos].DPorSP := 0;
        DPArray[DPNum] := bitpos;
        DPNum := DPnum +1;
      end;
      tp := tp shr 1;
    end;
  end;

  //对供应点设置距离为0，被供应的点加上本站。
  for j:=0 to SPNum -1 do
  begin
     MyGenes[SPArray[j]].Distance :=0;
     MyGenes[SPArray[j]].Supplypoint:= MyGenes[SPArray[j]].Supplypoint + inttostr(SPArray[j]);
  end;

  //对需求点求解他的供应点和距离；
  for j:=0 to DPNum -1 do
  begin
    tmpDistance :=1.0E+10;
    sp:=-1;
     for k := 0 to SPNum-1 do
     begin
       tmpX := power(Mygenes[SPArray[k]].x - Mygenes[DPArray[j]].x, 2) ;
       tmpY := power(Mygenes[SPArray[k]].y - Mygenes[DPArray[j]].y, 2) ;
       tmpDistance2 := sqrt(tmpX + tmpY);
       if (tmpDistance > tmpDistance2) then
       begin
         tmpDistance :=tmpDistance2;
         sp := SPArray[k];
       end;
     end;
     Mygenes[DPArray[j]].Distance := tmpDistance;
     if sp<> -1 then
       Mygenes[sp].Supplypoint:=Mygenes[sp].Supplypoint +',' + IntToStr(DPArray[j]);
  end;
  //计算目标值
  //计算建站成本；
   totalStationCost := 0;