工业控制之“什么叫RTO？”

读研究生时，过程控制领域经常涉及“APC”、“RTO”等字眼，导师也经常性提及，现在在工作中也开始提了，可能意识到先进控制的重要性了。

今天谈一下RTO在工业上的应用，曾经和一个博士生对“RTO涉及哪些算法”发生过激烈争吵，关于一些算法寻优名词，例如：”瞎子爬山法“、“进退法”、“成功失败法”等，这些算法其实是一种局部择优的贪心搜索算法，本质上是梯度下降法。神经网络中的权值调整就是通过梯度下降法。J目标函数（偏差最小），优化权值，典型案例求最值问题，梯度下降法先天不足就是只能在单调区间内寻优（单调函数），假设有多个单调区间，（负）梯度下降法就会被拉回原单调区间，这就是经常说的容易陷入局部最优解问题。权值w取值区间假设[0 +∞]，采用梯度下降法只能得到第一个谷点。如果应用到PID控制上，假设已知被控对象传递函数，经差分离散化得到差分函数，经过梯度下降法循环迭代求出一组权值（PID参数）得到第一个谷点，权值在迭代过程中始终没有越过[0 w5点]，也就是说限制在w3点上，典型案例人工神经网络优化PID参数。

关于梯度下降法原理可以参考：https://blog.csdn.net/weixin_43213884/article/details/122303703

在这里插入图片描述

不断逼近-0.625

遗传算法可以解决梯度下降法的局部最优解问题，举个完整例子[图片摘自智能控制导论_昆明理工大学]：

在这里插入图片描述

有几个自变量就有几个样本数据，把这几个样本数据合成一个染色体。

在这里插入图片描述

%Generic Algorithm for function f(x1,x2) optimum
clear all;
close all;

%Parameters
Size=80;   
G=100;     
CodeL=10;
 
umax=2.048;
umin=-2.048;
%初始化 种群 采样rand+round 产生0和1
E=round(rand(Size,2*CodeL));    %Initial Code

%Main Program
for k=1:1:G
time(k)=k;

for s=1:1:Size
m=E(s,:);
y1=0;y2=0;

%Uncoding  解码操作
m1=m(1:1:CodeL);
for i=1:1:CodeL
   y1=y1+m1(i)*2^(i-1);
end
x1=(umax-umin)*y1/1023+umin;
m2=m(CodeL+1:1:2*CodeL);
for i=1:1:CodeL
   y2=y2+m2(i)*2^(i-1);
end
x2=(umax-umin)*y2/1023+umin;
%优化目标是求函数的最大值，故可将个体的适应度直接取为对应的目标函数值
F(s)=100*(x1^2-x2)^2+(1-x1)^2;
end
%选个体适应度的倒数作为目标函数
Ji=1./F;     %（有一说一 这个真的没用）
%****** Step 1 : Evaluate BestJ ******
BestJ(k)=min(Ji);

fi=F;                          %适应度函数 
%'ascend' 表示升序（默认值），'descend' 表示降序。
[Oderfi,Indexfi]=sort(fi,'ascend' );     %升序排列
Bestfi=Oderfi(Size);           %Let Bestfi=max(fi )选出最优个体
BestS=E(Indexfi(Size),:);      %Let BestS=E(m), m is the Indexfi belong to max(fi)
bfi(k)=Bestfi;             %作图用的 寻优过程

%****** Step 2 : Select and Reproduct Operation******
   fi_sum=sum(fi);
   fi_Size=(Oderfi/fi_sum)*Size;%比例法进行 复制
   
   fi_S=floor(fi_Size);        %Selecting Bigger fi value
   
   kk=1;       
   for i=1:1:Size
      for j=1:1:fi_S(i)        %Select and Reproduce 
       TempE(kk,:)=E(Indexfi(i),:);  %%优秀的个体进行复制
         kk=kk+1;              %kk is used to reproduce
      end
   end
   
%************ Step 3 : Crossover Operation ************
pc=0.60;               %%交叉的概率
n=ceil(20*rand);       %%一点交叉 选择出一点交叉的位置
for i=1:2:(Size-1)
    temp=rand;
    if pc>temp                  %Crossover Condition
    for j=n:1:20
        TempE(i,j)=E(i+1,j);
        TempE(i+1,j)=E(i,j);
    end
    end
end
TempE(Size,:)=BestS;
E=TempE;
   
%************ Step 4: Mutation Operation **************
%pm=0.001;
%pm=0.001-[1:1:Size]*(0.001)/Size; %Bigger fi, smaller Pm
%pm=0.0;    %No mutation
pm=0.1;     %Big mutation   %%%变异

   for i=1:1:Size
      for j=1:1:2*CodeL
         temp=rand;
         if pm>temp                %Mutation Condition
            if TempE(i,j)==0
               TempE(i,j)=1;        %1变0 0变1
            else
               TempE(i,j)=0;
            end
        end
      end
   end
   
%Guarantee TempPop(30,:) is the code belong to the best individual(max(fi))
TempE(Size,:)=BestS;      %保证最优个体不丢失（变异交叉等会改变最优个体的基因）
E=TempE;                  %形成新的一代种群的基因
end
 
Max_Value=Bestfi
BestS
x1
x2
figure(1);
plot(time,BestJ); 
xlabel('Times');ylabel('Best J');
figure(2);
plot(time,bfi);
xlabel('times');ylabel('Best F');

在这里插入图片描述

关于遗传算法原理可以参考：https://blog.csdn.net/qq_30759585/article/details/112761360

遗传算法PID的解析和实现
https://zhuanlan.zhihu.com/p/139945063

function GA_PID()
%{
程序功能：
1、遗传算法PID的实现
2、传递函数：
G(s)=400/(s^2+50s)

date:2020.05.12
%}

    clear,clc
    close all

    global N  L   num w1 w2 w3 w4 vin t G
    
    %-----------------完成PID初始化---------------------------------------------
    
    w1=0.5;  %误差权重
    w2=0.4;  %超调量权重
    w3=0.9;  %上升时间
    w4=0.4;  %调节时间  %权重的大小可以影响对应的性能
     vin=1; %参考输入  
    
    KPRange=[0.001, 20];
    KIRange=[0.001 ,1];
    KDRange=[0.001,1 ] ;   %参数取值范围，不要从0开始，否则会陷入迭代困难！
    
    G=tf(400,[ 1 ,50, 0 ]) ; %被控对象
    time=0.1 ;%仿真时长
    t=linspace(0,time ,200);
    
    %----------------完成遗传算法初始化-----------------------------------------------

    num=3 ;%变量个数
    N = 50;                          % 种群上限
    ger = 10;                       % 迭代次数
    L = 10;                           % 单个变量的编码基因长度
    pc = 0.8;                        % 交叉概率
    pm = 0.1;                        % 变异概率
    pt=0.2;                         %自我复制概率

    dna1 = randi([0, 9], [N, L]);     % 基因
    dna2=randi([0, 9], [N, L]);     % 基因
    dna3=randi([0, 9], [N, L]);     % 基因

   
    
    dna=[dna1,dna2,dna3]; %合并3个变量的基因

    for epoch = 1: ger               % 进化代数为100
        %交叉
        [x1,x2]=jiaoCha(dna, pc);  
        % 变异操作
        x3=bianYi(dna, pm); 
        %自我复制
        x4=copy(dna, pt);
        dna = [dna; x1; x2; x3; x4];               % 合并新旧基因
        
        dna1=dna(:, 1:L);
        dna2=dna(:, L+1 : 2*L);
        dna3=dna(:, 2*L+1 : num*L);
        
        fi = zeros(N*5, 1);                % 初始化适应度，提速
%         fi = Fx(jieMa(dna1,xRange) , jieMa(dna2 ,yRange));    %         计算适应度，容易理解，求最大值
        M=length(dna1);
        for k=1:M
        fi(k)=-Fx(jieMa(dna1(k,:),KPRange) , jieMa(dna2(k,:) ,KIRange) , jieMa(dna3(k,:) ,KDRange)); %求最小值
        
        end
        
%         fi=Fx(jieMa(dna1,xRange) , jieMa(dna2 ,yRange)); %求最大值
        dna = [dna, fi];  %适应度放在最后一列  
        %自然选择
        dna=choice(dna, fi);
        dna = dna(:, 1 : num*L);  %去掉矩阵末尾的适应度，进入下一次迭代
        
      
    end
     x = jieMa( dna(:, 1:L) ,KPRange);            % 对最终种群解码 (dco(1)-1)
     y = jieMa( dna(:, L+1:2*L ) ,KIRange);
     z=  jieMa(dna( :, 2*L+1 : num*L) ,KDRange) ;

    KP=x(1)
    KI=y(1)
    KD=z(1)
%     Fxy=Fx(KP,KI,KD)
    y=draw( KP, KI, KD) ;
    figure
    plot(t ,y, 'linewidth' ,3) %绘制效果图
    title('Genetic algorithm PID')
    xlabel('time(s)')
    ylabel('Amplitude')
    grid on
    
    
end


%适应度函数/目标函数
function J=Fx(kp , ki , kd)
    global w1 w2 w3 w4 vin  
    [y, tr ,ts ,pos]=performance( kp, ki, kd ) ;
    n=length(y);
    J=0;
    for k=1:n
        e=vin-y(k) ; %误差
       
        J=J+w1*abs(e)+w2*pos+w3*tr ;      %误差超调量的积分  
        
    end
    J=J+w4*ts ;%调节时间
    
end

%绘制图形
function y=draw( kp, ki, kd)
    global t G
    Gpid=tf([kd, kp ,ki ], [1,0] ) ; %描述PID控制器的传递函数
    Gs=feedback(Gpid*G  ,1,  -1)  ;%描述闭环系统的负反馈传递函数，H=1
%     figure
    y=step(Gs, t) ; %单位阶跃响应
    

end

%{
计算PID性能指标
  tr:上升时间
  ts:调节时间
  pos:超调量
%}
function [y, tr ,ts ,pos]=performance(kp, ki, kd)
    
    global t  

	y=draw( kp, ki, kd);
    ymax=max(y);  %最大响应
    yend=y(length(y));  %响应终值
    pos=( ymax-yend )/yend;  %超调量
    
    n=1;
    while( y(n)<yend  )
        n=n+1;
    end
    tr=t(n) ;%上升时间
%     n=1;
%     while( y(n)<ymax )
%         n=n+1;
%     end
%     tp=t(n) ; %峰值时间
    n=length(t) ;
    while( y(n)>0.98*yend && y(n)<1.02*yend )
        n=n-1;
    end
    ts=t(n); %调节时间
    
    
       
end

%解码
function result=jieMa(gene , range)
     global L
     a=range(1) ;
     b=range(2) ;
    dco=zeros(L,1);
    for i=1 : L
        dco(i)=10^(L-i); %解码器
    end
    result=gene * dco / (dco(1)-1)/10 *(b-a)+a;
    
end



%选择
function result=choice(gene, fi)
    % 自然选择--排名法
    global L N num
    dna = flipud(sortrows(gene, L*num + 1));    % 对适应度进行降序排名，
    while size(dna, 1) > N                 % 自然选择,只保留50条基因
        d = randi(size(dna, 1));           % 排名法
        %当d非常大时，则容易被淘汰
        if rand < (d - 1) / size(dna, 1)  %根据概率大小，淘汰排名靠后的小概率基因 
            dna(d,:) = [];
            fi(d, :) = [];
        end
    end
    result=dna ; 
end

%变异
function result=bianYi(gene, pm)
    global N L num
    for i = 1: N                           % 变异操作
        if rand < pm
            gene(i,randi(num*L)) = randi([0, 9]);
        end
    end
    result=gene;

end

%复制
function result=copy(gene, pt)
    global N L num
    x= zeros(N, num*L);                % 初始化子代基因，提速用
   
    for i = 1: N                 % 交叉操作
        if rand < pt
           d = randi(N);                 
           x(i,:)=gene(d,:);% 随机抽选一个自我复制的个体
        end
    end
    result=x;
end

%交叉函数，一条dna，产生两条交叉后的dna
function [x1,x2]=jiaoCha(gene, pc)
    global N L num
    x1 = zeros(N, num*L);                % 初始化子代基因，提速用
    x2 = x1;                         % 同上
    for i = 1: N                 % 交叉操作
        if rand < pc
           d = randi(N);            
           m = gene(d,:);            % 随机抽选一个交叉的个体
           d = randi(num*L-1);          % 确定交叉断点
           x1(i,:) = [gene(i,1:d), m(d+1: num*L)];  % 新个体 1        
           x2(i,:) = [m(1:d), gene(i, d+1: num*L)];  % 新个体 2
        end
    end
end