采用实数编码的遗传算法进行PID参数整定的学习

遗传算法是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法。其计算过程主要分为六个步骤，详见百度百科。
本文主要是我在学习过程中对一段代码的理解的介绍。这段代码采用的是实数编码的方式（还有一种二进制编码方式），主要分为两部分：
一部分作为主函数用来循环迭代，进行选择，交叉，变异等步骤，其中部分属于固定操作；
另一部分是计算适应度函数的，这部分需要根据不同问题进行设计，学过运筹学的可以把适应度函数理解为目标函数，要整定的三个参数为变量，需要限制在一定范围内，参数整定就是在这范围内找到是目标函数取极值的变量大小。

其中主函数如下（注释已经详细表达了我的理解）：

clear all;
close all;
global rin yout timef

Size=30;  % 种群大小30个 可行解
CodeL=3;  % 需要辨识三个变量KP、KI、KD

%%%PID参数设置范围
MinX(1)=zeros(1);
MaxX(1)=20*ones(1);

MinX(2)=zeros(1);
MaxX(2)=1.0*ones(1);

MinX(3)=zeros(1);
MaxX(3)=1.0*ones(1);

%%%PID随机取值
Kpid(:,1)=MinX(1)+(MaxX(1)-MinX(1))*rand(Size,1);
Kpid(:,2)=MinX(2)+(MaxX(2)-MinX(2))*rand(Size,1);
Kpid(:,3)=MinX(3)+(MaxX(3)-MinX(3))*rand(Size,1);

G=100;    % 种群运行100次
BsJ=0;

%*************** Start Running ***************
for kg=1:1:G
    time(kg)=kg;%%%无用，记录循环次数

    %****** Step 1 : Evaluate BestJ ******
    for i=1:1:Size
        Kpidi=Kpid(i,:);  %%%提取单个个体（PID参数）

        [Kpidi,BsJ]=chap_f(Kpidi,BsJ);%%%输出Kpidi无用，BSJ为个体适应度

        BsJi(i)=BsJ;%%%记录一代30个体适应度
    end

    [OderJi,IndexJi]=sort(BsJi);%%%适应度排序
    BestJ(kg)=OderJi(1);%%%取最小适应度
    BJ=BestJ(kg);%%%无用
    Ji=BsJi+1e-10;    %%%避免分母为零

    fi=1./Ji;%%%取倒数，求极小变为求极大


    [Oderfi,Indexfi]=sort(fi);    %%%将极大目标值排序
    Bestfi=Oderfi(Size);          %%%取最大目标值，即最小适应度倒数
    BestS=Kpid(Indexfi(Size),:);  %%%最优个体对应的PID参数

%%%以下步骤是固定套路
    %****** Step 2 : 选择复制操作******
    fi_sum=sum(fi);
    fi_Size=(Oderfi/fi_sum)*Size;%%%归一化目标函数值，再乘种群

    fi_S=floor(fi_Size);         %%%floor化整数
    r=Size-sum(fi_S);

    Rest=fi_Size-fi_S;
    [RestValue,Index]=sort(Rest);

    for i=Size:-1:Size-r+1
        fi_S(Index(i))=fi_S(Index(i))+1;     % Adding rest to equal Size
    end

    k=1;
    for i=Size:-1:1       % Select the Sizeth and Reproduce firstly
        for j=1:1:fi_S(i)
            TempE(k,:)=Kpid(Indexfi(i),:);      % Select and Reproduce
            k=k+1;                            % k is used to reproduce
        end
    end

    %************ Step 3 : 交叉操作************
    Pc=0.90;%%%交叉概率
    for i=1:2:(Size-1)
        temp=rand;
        if Pc>temp                      %%%相邻两个以一个随机数比例分配变化
            alfa=rand;
            TempE(i,:)=alfa*Kpid(i+1,:)+(1-alfa)*Kpid(i,:);
            TempE(i+1,:)=alfa*Kpid(i,:)+(1-alfa)*Kpid(i+1,:);
        end
    end
    TempE(Size,:)=BestS;%%%把最好的个体保留
    Kpid=TempE;

    %************ Step 4: 变异操作 **************
    Pm=0.10-[1:1:Size]*(0.01)/Size;       %%%最大目标值（最小适应度）对应最小的变异概率，向量
    Pm_rand=rand(Size,CodeL);
    Mean=(MaxX + MinX)/2; %%%取范围中值
    Dif=(MaxX-MinX);%%%取范围大小

    for i=1:1:Size
        for j=1:1:CodeL
            if Pm(i)>Pm_rand(i,j)        %%%变异，范围中值加上总范围值乘(rand-0.5)
                TempE(i,j)=Mean(j)+Dif(j)*(rand-0.5);
            end
        end
    end
   
    TempE(Size,:)=BestS;
    Kpid=TempE;
end
Bestfi
BestS
Best_J=BestJ(G)
figure(1);
plot(time,BestJ);
xlabel('Times');ylabel('Best J');
figure(2);
plot(timef,yout,'b');
xlabel('Time(s)');ylabel('rin,yout');

求适应度如下

function [Kpidi,BsJ]=chap_f(Kpidi,BsJ)
global rin yout timef

ts=0.001;
sys=tf(400,[1,50,0]);  %%%定传递函数，并离散
dsys=c2d(sys,ts,'z');
[num,den]=tfdata(dsys,'v');

rin=1.0;%%%设定值
u_1=0.0;u_2=0.0;
y_1=0.0;y_2=0.0;
x=[0,0,0]';%%%PID数字化对应需要用到的偏差的计算值
B=0;%%%适应度，即需要极小化的目标函数
error_1=0;
tu=1;%%%上升时间，目标函数中的东西
s=0;
P=100;


for k=1:1:P
    timef(k)=k*ts;
    r(k)=rin;
    %%%计算输入U
    u(k)=Kpidi(1)*x(1)+Kpidi(2)*x(2)+Kpidi(3)*x(3);
    %%%限制U
    if u(k)>=10
        u(k)=10;
    end
    if u(k)<=-10
        u(k)=-10;
    end
    %%%计算输出Y
    yout(k)=-den(2)*y_1-den(3)*y_2+num(2)*u_1+num(3)*u_2;
    %%%计算偏差
    error(k)=r(k)-yout(k);
    %------------ Return of PID parameters -------------
    u_2=u_1;u_1=u(k);
    y_2=y_1;y_1=yout(k);


    x(1)=error(k);                % Calculating P
    x(2)=(error(k)-error_1)/ts;   % Calculating D
    x(3)=x(3)+error(k)*ts;        % Calculating I
    error_1=error(k);

    if s==0
        if yout(k)>0.95&&yout(k)<1.05
            tu=timef(k); % 上升时间
            s=1;
        end
    end
end

for i=1:1:P
    Ji(i)=0.999*abs(error(i))+0.01*u(i)^2*0.1;%%%限制偏差和U
    B=B+Ji(i);
    if i>1
        erry(i)=yout(i)-yout(i-1);
        if erry(i)<0
            B=B+100*abs(erry(i));%%%限制超调
        end
    end
end
BsJ=B+0.2*tu*10;%%%输出个体适应度
end

运行结果如下，但我发现每次运行结果都有所不同，但大致都差不多，我理解为可能是其中的rand函数导致的吧。