MATLAB遗传算法PID大作业.doc

Four short words sum up what has lifted most successful individuals above the crowd: a little bit more.-author-dateMATLAB遗传算法PID大作业关于加强课程考试试题、试卷、标准答案规范要求遗传算法在调节控制系统参数中的应用【摘要】自动化控制系统多采用PID控制器来调节系统稳定性和动态性，PID的Kp,Ki,Kd参数需要合理选择方能达到目标。遗传算法是一种模拟生物进化寻求最优解的有效算法，本文通过利用GAbx工具箱实现对控制电机的PID进行参数优化，利用matlab的仿真功能可以观察控制效果。1. 直流伺服电机模型1.1物理模型图1 直流伺服电机的物理模型-电枢输入电压（） -电枢电阻（） -电枢电感（H）-感应电动势（） -电机电磁转矩（N） J-转动惯量（） B-粘性阻尼系数（） -流过电枢的电流（A） -电机输出的转角（）1.2传递函数利用基尔霍夫定律和牛顿第二定律得出电机基本方程并进行拉布拉斯变换式中：为电机的转动常数（）；为感应电动势常数（）图2 直流伺服电机模型方框图消去中间变量得系统的开环传递函数：系统参数如下：2. PID校正图3 PID校正 Kp,Ki,Kd为比例，积分，微分系数令Kp=15、Ki=0.8 、Kd=0.6M文件：J=3.23E-6;B=3.51E-6;Ra=4;La=2.75E-6;Kt=0.03;num= Kt;den=(J*La) (J*Ra)+(La*B) (B*Ra)+Kt*Kt) 0;t=0:0.001:0.2;step(num,den,t);Kp=15;Ki=0.8;Kd=0.6;numcf=Kd Kp Ki;dencf=1 0;numf=conv(numcf,num);denf=conv(dencf,den);numc,denc=cloop(numf,denf);t=0:0.001:0.04;step(numc,denc,t);matlab进行仿真，我们可以看出不恰当的PID参数并不能使系统达到控制系统的要求，因此需要对PID参数进行优化。图4 系统阶跃响应（Kp=15，Ki=0.8，Kd=0.6）3. 遗传算法3.1 遗传算法和工具箱简介遗传算法（GA)是基于自然选择和基因遗传学原理的优化搜索方法。它借鉴了达尔文的进化论和孟德尔的遗传学说。其本质上是一种高效、并行、全局搜索的方法，它能在搜索过程中自动获取和积累有关搜索空间的知识，并自适应地控制搜索过程以求得最优解。遗传算法操作使用适者生存的原则，在潜在的解决方案种群中逐次产生一个近似最优的方案。在遗传算法的每一代中，根据个体在问题域中的适应度值和从自然遗传学中借鉴来的再造方法进行个体选择，产生一个新的近似解。在这个过程导致种群中个体的进化，得到的新个体比原个体更能适应环境，就像自然界中的改造一样。表3.1遗传学和遗传算法中基本用语对照表遗传学遗传算法染色体（Chromosome)解的编码（算法的操作对象）基因（Gene）解中每一分量等位基因（Allele）特性值基因座（Locus）二进制串中位置基因型（Genptype）结构表现型（Phenotype）参数集、候选解个体（Individual）解适者生存在算法停止时，最优目标值的解有最大可能被留住适应性（Fitness）适应度函数值群体（Population）选定的一组解复制（Reproduction）根据适应度函数值选取的一组解交配（Crossover）通过交配产生一组新解的过程变异（Mutation）编码的某一个分量发生变化的过程英国谢菲尔德大学开发的遗传算法工具箱把参数，选择，交叉，变异等过程封装成函数进行操作，其基本搜索过程不变。表3.2 遗传算法工具箱常用函数创建种群crtbase创建基向量crtbp创建任意离散随即种群crtrp创建实值初始种群（bs2rv）适应度计算ranking常用的基于秩的适应度计算scaling比率适应度计算选择函数reins一致随机和基于适应度的重插入rws轮盘选择select高级选择例程sus随机遍历采样变异算子mut离散变异mutate高级变异函数mutbga实值变异交叉算子recdis离散重组recint中间重组reclin线性重组recmut具有变异特征的线性重组recombine高级重组算子xovdp两点交叉算子xovdprs减少代理的两点交叉xovmp通常多点交叉xovsh洗牌交叉xovshrs减少代理的洗牌交叉xovsp单点交叉xovsprs减少代理的单点交叉子种群的支持migrate在子种群间交换个体实用函数bs2rv二进制串到实值的转换rep矩阵的复制3.2 利用遗传算法优化过程1）根据遗传算法优缺点这里选择遗传代数为100，种群大小为30，变量维数为3其中10Kp20,0Ki1, 0Kd1，要求精度0.0001所以二进制串编码长度为17故Kp精度为(20-10)/( -1) Ki,Kd精度为（1-0）/(-1),代沟为0.9，交叉概率为0.6，变异概率为0.01。2）适应度函数的设计PID优化设计的目的是使系统某些性能指标最优，然而，单纯的误差性能指标很难同时满足系统对快速性、稳定性和鲁棒性的要求，因此在适应度函数中引入超调量、上升时间和累计绝对误差指标项。设=(|yp-yref|/yref)×100%为系统的超调量，yp和yref分别为输出峰值和输入参考值;t*r为上升时间，将其定义为输出从0第1次达到0.95yref的时间;ek=yk-yref为采样时刻k的输出误差;wj(j=1，2，3)为权重系数。则多目标适应度函数为: (5)通过对权重系数的调整，可以改变系统对快速性和稳定性的要求。如系统要求较小的超调，可以适当增大w1;若系统要求快速的动态响应，则可以适当增大w2。计算机控制是一种采样控制，它只能根据采样时刻的偏差值计算控制量。因此连续PID控制算法不能直接使用，需要采用离散化方法。在计算机PID控制中，使用的是数字PID控制器。图3.2 增量式PID控制系统增量式PID控制是数字控制器的输出只是控制器的增量u(k)。当执行机构需要的是控制量的增量时，应采用增量式PID控制，采样时刻点kT（T为采样周期）代表连续时间t。做如下近似变换来离散化：可得离散的PID表达式： 增量式PID控制算法： 增量式PID控制算法不需多次累加，控制增量u(k)仅与最近3次得采样有关，所以误动作时的影响相对较小图3.3 基于GA算法参数调整的PID控制系统采样时间取1ms，为获取满意的过渡过程动态特性，采用误差绝对值时间积分性能指标作为参数选择的目标函数。为防止控制量过大，在目标函数中加入控制输入的平方项。选用下式作为参数选取的最优指标 为系统偏差，为PID控制器输出，为上升时间。 为避免超调了，采用惩罚功能，一旦产生超调，将超调量作为最优指标的一项，此时最优指标为 J值越小，对应的PID参数越优良。但是，由于GA算法用于求解最大值问题，所以需要对J做一个简单变形，从而得到算法的适应度函数：3）依据遗传算法的步骤编写程序并仿真图3.4 BestKpKiKd优化值时的阶跃响应图3.5 J随进化代数的变化曲线图3.6 优化得到BestKpKiKd及此时的LeastJ4. 程序代码1.GA.m文件%GA(Generic Algorithm) Program to optimize Parameters of PIDclcclear all;close all;             %清除命令，变量，关闭窗口global rin yout timef  %输入，输出，离散时间%*Initialization*%MAXGEN=100;            %遗传代数NIND=30;               %种群大小NVAR=3;                %变量维数PRECI=17;              %变量精度GGAP=0.9;              %代沟(Generation gap)Pc=0.6;                %交叉概率Pm=0.01;               %变异概率FieldD=PRECI,PRECI,PRECI;10,0,0;20,1,1; 1,1,1;0,0,0;1,1,1;1,1,1;   %区域描述器       Chrom=crtbp(NIND,NVAR*PRECI);                                       %初始种群 for i=1:1:MAXGEN    gen(i)=i;                         %遗传代数   %* Step 1 : Evaluate LeastJ *      KpKiKdlsj=bs2rv(Chrom,FieldD);     %各代种群由二进制串转化为实值  for j=1:1:NIND    J=0;    KpKiKd=KpKiKdlsj(j,:);                    KpKiKd,J=pidf(KpKiKd,J);        %计算本代种群中各个个体的目标函数值    Lsj(j,1)=J;   end  OderLsj,IndexLsj=sort(Lsj);       %本代种群中各个个体的目标函数值升序排列  genlstj(i)=OderLsj(1);              %升序排列中第一个最小，为本代最优  BestKpKiKd=KpKiKdlsj(IndexLsj(1),:);%BestKpKiKd    Lsj=Lsj+1e-10;                      %Avoiding deviding zero 1*10(-10)   Msj=1./Lsj;                         %设置适应度函数       OderMsj,IndexMsj=sort(Msj);       %Arranging fi small to bigger适应度值升序排列  genmost=OderMsj(NIND)   ;           % Let Bestfi=max(fi)升序中最后一个最大，为本代最优  moststring=Chrom(IndexMsj(NIND),:);    %本代中对应最优的PID参数 %* Step 2 : Select and Reproduct Operation*  FitnV = ranking(Lsj);               %分配适应度值(Assign fitness values)             SelCh=select('sus',Chrom,FitnV,GGAP);  %选择 %* Step 3 : Crossover Operation *  SelCh=recombin('xovsp', SelCh,Pc);     %交叉   %* Step 4: Mutation Operation *  SelCh=mut(SelCh,Pm);                   %变异  SelCh(NIND,:)=moststring;           %保证交叉后的新种群Temp2包含原种群中的最优个体%*Step 5：Reinsert Operation*  ObjV=FitnV;     KpKiKdlsj=bs2rv(SelCh,FieldD);      %计算ObjVSel  for j=1:1:NIND    J=0;    KpKiKd=KpKiKdlsj(j,:);                    KpKiKd,J=pidf(KpKiKd,J);            Lsj(j,1)=J;   end  Lsj=Lsj+1e-10;                         Msj=1./Lsj;   ObjVSel=Msj;    Chrom ObjV = reins(Chrom, SelCh, 1, 1, ObjV,ObjVSel);    %重插入 end  LeastJ=genlstj(MAXGEN)              %最优目标函数值（性能指标） BestKpKiKd                          %最优PID参数（实值） figure(1); plot(gen,genlstj); xlabel('gen(s)');ylabel('genlstj'); figure(2); plot(timef,rin,'r',timef,yout,'b'); xlabel('Time(s)');ylabel('rin,yout');2.pidf.m函数文件function KpKiKd,J=pidf(KpKiKd,J)global rin yout timef tz=0.001;sys=tf(0.03,8.8825e-12,1.2920e-05,9.1404e-04,0); zsys=c2d(sys,tz,'z');num,den=tfdata(zsys,'v'); rin=1.0;u_1=0.0;u_2=0.0;u_3=0.0;y_1=0.0;y_2=0.0;y_3=0.0;K=0,0,0'error_1=0;Tup=1;m=0;final=100; for t=1:1:final timef(t)=t*tz; r(t)=rin; u(t)=KpKiKd(1)*K(1)+KpKiKd(2)*K(2)+KpKiKd(3)*K(3); if u(t)>=10 u(t)=10; end if u(t)<=-10 u(t)=-10; end yout(t)=-den(2)*y_1-den(3)*y_2-den(4)*y_3+num(2)*u_1+num(3)*u_2+num(4)*u_3; error(t)=r(t)-yout(t);%- Return of PID parameters - u_3=u_2;u_2=u_1;u_1=u(t); y_3=y_2;y_2=y_1;y_1=yout(t); K(1)=error(t); % Calculating P K(2)=(error(t)-error_1)/tz; % Calculating D K(3)=K(3)+error(t)*tz; % Calculating I error_1=error(t);if m=0 if yout(t)>0.95&yout(t)<1.05 Tup=timef(t); m=1; end endend for t=1:1:final J=J+0.999*abs(error(t)+0.001*u(t)2; if t>1 erry(t)=yout(t)-yout(t-1); if erry(t)<0 J=J+100*abs(erry(t); end endendJ=J+2*Tup;-