velocity.~pas

这是一个大学里有关粒子群算法的程序

unit velocity;

interface

const nDim = 3;           // 维度
      popSize = 20;       // 粒子群大小
      prev = 100;         // 随机数精度
      presion = 1e-20;    // 适应值退出精度要求
type

  Tvelocity = class(TObject)
    protected
      c1,c2,w,fpBest,fBestVal : double;
      curPoint,Vmax : Array [0..nDim-1] of double;
      pRand : Array [0..nDim-1,0..1] of double;
      VelPostion,curVelo,curPostion : Array[0..popSize-1,0..nDim-1] of double;
      FValue,cValue : Array [0..popSize-1] of double;
      pBest : Array[0..nDim-1] of double;
      pBestIdx,gBestIdx,iteration,maxIter : Integer;
      oldBestVal: INT64;


    public
      FgBest: Array[0..nDim-1] of double;

    private
      function func1(x,y,z:double):double;   //方程一
      function func2(c,d,e:double):double;   //方程二
      function func3(m,n,q:double):double;   //方程三

    public
      Constructor Create;

      function mainFun:string;
      procedure setRand(nDim:Integer; min,max:double);
      procedure genRandData;
      procedure genVmax;
      function genFValue:double;

  end;

implementation

Uses math,SysUtils,Dialogs;


Constructor Tvelocity.Create;
begin
  oldBestVal := 0;            // 相关初始化工作
  iteration := 0;            // 迭代序号
  maxIter := 200;            // 最大迭代次数

  c1 := 0.5;
  c2 := 0.5;
  w := 0.7;

  setRand(0,3.0,5.0);         // 设置粒子群各维度的上下界
  setRand(1,2.0,4.0);
  setRand(2,0.5,2.0);
end;

procedure TVelocity.genRandData;
var
  i : Integer;
begin
  for i := 0 to nDim-1 do
    curPoint[i] := random(popSize*prev)*(pRand[i][1]-pRand[i][0])
        /popSize/prev+pRand[i][0];
end;

procedure TVelocity.genVmax;
var
  i : Integer;
  min,max : double;
begin
  for i := 0 to nDim-1 do
  begin
    max := pRand[i][1];
    min := pRand[i][0];
    Vmax[i] := (max-min)/2;
  end;
end;

procedure TVelocity.setRand(nDim:Integer; min,max:double);
begin
  pRand[nDim][0] := min;
  pRand[nDim][1] := max;
end;

function TVelocity.func1(x,y,z:double):double;   //方程一
begin
   result:=power(x,y)+power(y,x)-5*x*y*z-85;
end;

function TVelocity.func2(c,d,e:double):double;   //方程二
begin
   result:=power(c,3.0)-power(d,e)-power(e,d)-60;
end;

function TVelocity.func3(m,n,q:double):double;   //方程三
begin
   result:=power(m,q)+power(q,m)-n-2;;
end;

function TVelocity.genFValue : double;
begin
  // 增加计算适应度值
  result := func1(curPoint[0],curPoint[1],curPoint[2])
          *func1(curPoint[0],curPoint[1],curPoint[2])
          +func2(curPoint[0],curPoint[1],curPoint[2])
          *func2(curPoint[0],curPoint[1],curPoint[2])
          +func3(curPoint[0],curPoint[1],curPoint[2])
          *func3(curPoint[0],curPoint[1],curPoint[2]);
end;

function Tvelocity.mainFun:string;
var
  isStop,i,j : Integer;
  minValue,tempValue,gBestValue : double;
  R1,R2 : Array [0..nDim-1] of double;
begin
  isStop := 0;
  Randomize;
  for i := 0 to popSize-1 do  // 生成粒子种群
  begin
    genRandData;              // 生成一个粒子
    VelPostion[i][0] := curPoint[0];        // 保存粒子
    VelPostion[i][1] := curPoint[1];
    VelPostion[i][2] := curPoint[2];
    FValue[i] := genFValue;         // 计算CurPoint粒子的适应值并保存
  end;

  genVmax;      // 生成粒子种群中各点的最大速度
  for i := 0 to popSize-1 do  // 生成粒子群中各点的速度
    for j := 0 to nDim - 1 do
      curVelo[i][j] := random(popSize*prev)*Vmax[j]/popSize/prev;

  pBestIdx := 0;
  minValue := FValue[0];
  for i := 1 to popSize-1 do       //  查找第一个群内的最优粒子
  begin
    if minValue > FValue[i] then
    begin
      minValue := FValue[i];
      pBestIdx := i;
    end;
  end;                            // minValue 为当前最优粒子的适应值

  pBest[0] := VelPostion[pBestIdx][0];    // 保存当前群最优粒子
  pBest[1] := VelPostion[pBestIdx][1];
  pBest[2] := VelPostion[pBestIdx][2];

  gBestValue := minValue;
  FgBest[0] := VelPostion[pBestIdx][0];    // 保存全局最优粒子
  FgBest[1] := VelPostion[pBestIdx][1];
  FgBest[2] := VelPostion[pBestIdx][2];
  gBestIdx := pBestIdx;

  while (isStop = 0) and (iteration < MaxIter) do   // 如相同的值超过一定次数则退出
  begin
    inc(iteration);                 // 调整迭代次数

    for i := 0 to nDim-1 do         // 生成随机数R1、R2
    begin
      R1[i] := random(popSize*prev)/popSize/prev;
      R2[i] := random(popSize*prev)/popSize/prev;
    end;

    for i := 0 to popSize-1 do            // 生成新群，并处理
    begin
      for j := 0 to nDim - 1 do
      begin
        curVelo[i][j] := w*curVelo[i][j]+c1*R1[j]*(pBest[j]-VelPostion[i][j])
              +c2*R2[j]*(FgBest[j]-VelPostion[i][j]);      // 调整速度
        curPoint[j] := VelPostion[i][j]+curVelo[i][j];    // 生成新的粒子
        while curPoint[j] >= pRand[j][1] do         // 上界越界后调整
          curPoint[j] := curPoint[j] - abs(CurVelo[i][j]);
        while curPoint[j] <= pRand[j][0] do         // 下界越界后调整
          curPoint[j] := curPoint[j] + abs(CurVelo[i][j]);
      end;
      tempValue := genFValue;
      if tempValue < FValue[i] then           // 新粒子优于原位粒子时，修改粒子
      begin
        VelPostion[i][0] := curPoint[0];      // 修改粒子
        VelPostion[i][1] := curPoint[1];
        VelPostion[i][2] := curPoint[2];
        FValue[i] := tempValue;               // 修改适应值
      end;
    end;
    pBestIdx := 0;
    minValue := FValue[0];
    for i := 1 to popSize-1 do       //  查找当前群内的最优粒子
    begin
      if minValue > FValue[i] then
      begin
        minValue := FValue[i];
        pBestIdx := i;
      end;
    end;
    pBest[0] := VelPostion[pBestIdx][0];    // 保存当前最优粒子
    pBest[1] := VelPostion[pBestIdx][1];
    pBest[2] := VelPostion[pBestIdx][2];

    if gBestValue > minValue then   // 全局适应值大于当前最优的，修改全局适应值
    begin
      FgBest[0] := VelPostion[pBestIdx][0];
      FgBest[1] := VelPostion[pBestIdx][1];
      FgBest[2] := VelPostion[pBestIdx][2];
      gBestValue := minValue;
      gBestIdx := pBestIdx;
    end;

    if gBestValue < presion then          // 精度达到要求，即退出
      isStop := 1

  end;
  result := 'x1 := '+FloatToStr(FgBest[0])+'    x2 := '+FloatToStr(FgBest[1])
          +'    x3 := '+FloatToStr(FgBest[2])+'        FValue = '+FloatToStr(gBestValue)
          +'    迭代次数：'+IntToStr(iteration);
end;



end.