MWORKS电力电子与电机系统建模与仿真 课件 第六章 基于MWORKS的最优化算法与群优化算法

第六章基于MWORKS的最优化算法与群优化算法目录6.1一维搜索算法6.2线性优化算法6.3蚁群优化算法6.4PSO6.1一维搜索算法最优化理论最优化问题通常可以表示为数学规划问题，至今已出现线性规划、整数规划、非线性规划、几何规划、动态规划、随机规划、网络流等许多分支，最优化理论和算法在实际应用中正在发挥越来越大的作用。一维搜索算法一维搜索算法也称为线性搜索法，是多变量求解算法的一部分，因此受到普遍重视。6.1一维搜索算法数学规划是指对含有几个变量的目标函数求极值，而这些变量也可能受到某些条件（等式方程或不等式方程）的限制，其一般数学表达式为：决策变量n维向量x=(x1,x2,…,xn)T。目标函数又称成本函数，为f

(x)。可行集集合Ω是n维实数空间Rn的一个子集，称为约束集或可行集。该优化问题的含义是，寻找合适的x，使得函数

f达到最小。6.1一维搜索算法上述问题是有约束优化问题的一般形式，约束变量只能从约束集Ω中取值。如果Ω∈Rn，则该问题是一个无约束优化问题。如果是一个有约束优化问题，则约束集可表示为Ω={x:h(x)=0,g(x)≤0}，其中，h(x)为等式约束条件，g(x)为不等式约束条件。讨论优化问题的最优性条件需要给出两类极小点的定义。6.1一维搜索算法【定义6.1】存在一个n元实值函数f:Rn→R，定义域。对于定义域Ω中的一个点x*，如果存在ε>0，对于所有满足||x-x*||<ε，||x-x*||<ε的向量x，不等式f

(x)≥f

(x*)都成立，则称x*是函数f的定义域Ω中的一个局部极小点；如果对于所有x∈Ω\{x*}，不等式f

(x)≥f

(x*)都成立，则称x*是函数f的定义域Ω中的一个全局极小点。6.1一维搜索算法如果将定义6.1中的f

(x)≥f

(x*)替换为f

(x)>f

(x*)，那么局部极小点和全局极小点对应成为严格局部极小点与严格全局极小点。严格地说，只有得到优化问题的全局极小点，才说明优化问题得到了解决。但实际上，复杂的优化问题很难获得全局极小点。因此在实际应用中，一般将局部极小点作为优化问题的解值。6.1一维搜索算法黄金分制法可用于求解一元单值函数f:R→R在闭区间[a0,b0]上的极小点。用黄金分制法、斐波那契数列法和二分法求解该问题的唯一前提是目标函数f在区间[a0,b0]上是单峰的，即存在唯一的局部极小点。6.1一维搜索算法思路挑选区间[a0,b0]中的点，计算对应的目标函数值，通过比较不断缩小极小点所在的区间。如何选择合适的点？利用尽可能少的计算次数找出函数的极小点，即不断压缩极小点所在的区间，直到达到足够的精度水平。黄金分制法黄金分割点是指将一条线段分割为两部分，使得其中一部分与全长的比等于另一部分与这部分的比，其比值约为0.618。利用黄金分制法开展区间压缩，将区间的压缩比例因子设定为ρ=0.618。经过N步压缩，极小点所在区间的长度将被压缩到初始长度的(0.618)N，这称为总压缩比。6.1一维搜索算法STEP01确定试探点在搜索区间[a0,b0]上确定两个试探点：

a1=a0+(1−ρ)(b0−a0)

b1=a0+ρ(b0−a0)STEP02计算函数值分别计算这两个试探点的函数值f

(a1)和f

(b1)STEP03比较并缩小区间若f

(a1)<f

(b1)，则区间[b1,b0]不可能存在极小点，区间范围缩小为[a0,b1]，即更新搜索范围边界值b0=b1；否则，区间[a0,a1]不可能存在极小点，区间范围缩小为[a1,b0]，即更新搜索范围边界值a0=a1STEP04精度判断若新的区间[a0,b0]小于搜索精度ε，则搜索停止；否则，回到步骤（1）继续搜索6.1一维搜索算法▌例6.1利用黄金分割法确定函数f

(x)=x4−14x3+60x2−70x在区间[0,2]中的极小点，要求将极小点所在的区间范围压缩到0.01之内。使用Julia语言在Syslab中编程，代码如下：clear()g(x)=x^4-14*x^3+60*x^2-70*x#定义目标函数a=0#定义搜索区间b=2h=b-aepsilon=0.01 #定义收敛要求rho=0.618 #定义黄金分割点a_1=a+(1-rho)*h #左搜索点b_1=a+rho*h #右搜索点6.1一维搜索算法n=1#迭代次数统计#迭代搜索whileh>epsilonglobala,b,a_1,b_1,h,f_1,f_2,n#计算左、右搜索点的函数值

f_1=g(a_1)f_2=g(b_1)iff_1<f_2#更新搜索区间

b=b_1h=b-ab_1=a_1a_1=a+(1-rho)*h

elsea=a_1h=b-aa_1=b_1b_1=a+rho*hendn=n+1#统计搜索次数endprint("x_min=$a_1,f_min=$f_1,InterationTime=$n")6.1一维搜索算法输出：x_min=0.7802345887428408,f_min=-24.369314188973373,InterationTime=13常规打印输出时，数据结果显示的小数位较多，可以通过sprintf或round语句实现指定小数位长度的数据格式。具体代码如下：x_min=@sprintf"%.4f"a_1f_min=round(f_1,digits=4)print("x_min=$x_min,f_min=$f_min,InterationTime=$n")输出：x_min=0.7802,f_min=-24.3693,InterationTime=136.1一维搜索算法区别：黄金分制法的求解过程中，压缩比ρ始终保持不变。斐波那契数列法允许参数ρ不断变化，获得比黄金分割法更高的压缩比。递推关系斐波那契数列是线性递推数列，在优化算法中设定初始项：F₀=0,F₁=1递推公式：Fₖ₊₁=Fₖ+Fₖ₋₁数列前15项F1F2F3F4F5F6F7F812358132134F9F10F11F12F13F14F155589144233377610987斐波那契数列法6.1一维搜索算法斐波那契数列法的总压缩比为：斐波那契数列法的总压缩比为：即最后一次迭代的参数为0.5。此时两个搜索点重合，无法对区间进行压缩。因此将最后一次迭代的参数设定为:修正后斐波那契数列法的总压缩比为：6.1一维搜索算法STEP01初始化给定初始搜索区间[a0,b0]和搜索精度ε，确定斐波那契数列的长度n，使得STEP02确定压缩比STEP03确定两个试探点a1=a0+ρ(b0−a0)b1=a0+(1−ρ)(b0−a0)并分别计算这两个试探点的函数值f

(a1)和f

(b1)ρ=1−F(n)/F(n+1)6.1一维搜索算法STEP04STEP05精度判断若新的区间[a0,b0]小于搜索精度ε，则搜索停止；否则，回到步骤（1）继续搜索比较并缩小区间若f

(a1)<f

(b1)，则区间[b1,b0]不可能存在极小点，区间范围缩小为[a0,b1]，即更新搜索范围边界值b0=b1；否则，区间[a0,a1]不可能存在极小点，区间范围缩小为[a1,b0]，即更新搜索范围边界值a0=a16.1一维搜索算法▌例6.2利用斐波那契数列法确定函数f

(x)=x4−14x3+60x2−70x在区间[0,2]中的极小点，要求将极小点所在区间范围压缩到0.01之内。使用Julia语言在Syslab中编程，代码如下：#Fibonaccimethodclear()g(x)=x^4-14*x^3+60*x^2-70*x#定义目标函数a=0 #定义搜索区间b=2epsilon=0.01 #定义收敛要求delta=0.01 #随机小量range=round((b-a)/epsilon,RoundUp)#向上取整Fi=zeros(50)Fi[1]=1Fi[2]=1n=1whileFi[n+1]<=range#确定迭代次数

globaln,FiFi[n+2]=Fi[n+1]+Fi[n]n=n+1endtime=n-1rou=1-Fi[n]/Fi[n+1]a_1=a+rou*(b-a) #左搜索点a_2=a+(1-rou)*(b-a) #右搜索点#迭代搜索6.1一维搜索算法whilen>1globala,b,a_1,a_2,f_1,f_2,n,rouf_1=g(a_1)#左搜索点的函数值

f_2=g(a_2)#右搜索点的函数值

n=n-1rou=1-Fi[n]/Fi[n+1]iff_1<f_2#更新搜索区间

b=a_2a_2=a_1f_2=g(a_2)ifn>2a_1=a+rou*(b-a)f_1=g(a_1)elsea_1=a_2-deltaf_1=g(a_1)endelsea=a_1a_1=a_2f_1=g(a_1)ifn>2a_2=a+(1-rou)*(b-a)f_2=g(a_2)elsea_2=a_1+deltaf_2=g(a_2)endendendx_min=round(a_1,digits=4)f_min=round(f_1,digits=4)print("x_min=$x_min,f_min=$f_min,InterationTime=$time")输出：x_min=0.7711,f_min=-24.3666,InterationTime=116.1一维搜索算法区别：二分法使用的是目标函数的导数f'

(x)，而不是黄金分割法和斐波那契数列法使用的函数值；在每次迭代中，区间压缩比为1/2，总压缩比为(1/2)N，其压缩比低于黄金分割法和斐波那契数列法的压缩比。基本思想二分法是一种最简单的试探法，其基本思想是，通过计算函数的导数值来缩短搜索区间。二分法本质上是求解目标函数的一阶导数f‘

(α)=0，因此要求函数f是一阶连续可微的。二分法6.1一维搜索算法01.确定初始中点给定初始搜索区间[a0,b0]，确定初始区间的中点：c=(a0+b0)/202.计算导数计算函数f在中点c处的一阶导数f'

(c)03.比较并缩小区间若f'

(c)>0，则说明极小点位于c点左侧，极小点所在区间可压缩为[a0,c]；若f'

(c)<0，则说明极小点位于c点右侧，极小点所在区间可压缩为[c,b0]04.精度判断若f'

(c)=0，则说明c点就是函数f的极小点，搜索停止；否则，回到步骤（1）继续搜。6.1一维搜索算法▌例6.3利用二分法确定函数f

(x)=x4−14x3+60x2−70x在区间[0,2]中的极小点，要求将极小点所在区间范围压缩到0.01之内。使用Julia语言在Syslab中编程，目标函数的一阶导数通过调用TySymbolicMath函数库中的derivative()函数生成，导数值的获取需要调用substitute()函数，代码如下：clear()usingTySymbolicMath@variablesxf=x^4-14*x^3+60*x^2-70*x#定义目标函数df=derivative(f,x)#定义目标函数的一阶导数a=0#定义搜索区间b=2epsilon=0.01#定义收敛要求6.1一维搜索算法n=1#迭代次数统计B=1#迭代搜索while(b-a)>epsilonglobala,b,c,n,xmin,B,y,dfc=(a+b)/2#计算区间中点

y=value(substitute(df,Dict(x=>c)))ify<0a=celseify==0xmin=cB=0breakelseb=cendendn=n+1endifB==1xmin=(a+b)/2endtime=n-1x_min=round(xmin,digits=4)fmin=value(substitute(f,Dict(x=>x_min)))#获取函数极小值f_min=round(fmin,digits=4)print("x_min=$x_min,f_min=$f_min,InterationTime=$time")输出：x_min=0.7773,f_min=-24.3692,InterationTime=86.1一维搜索算法设目标函数f

(x)是二阶连续可微的，将函数f

(x)在迭代点x(k)处泰勒展开，并取二阶近似：求解函数f

(x)的极小点可近似求解g(x)的极小点，

g(x)的极小点应满足一阶必要条件:令x=x(k+1)，可得基本思想若目标函数连续二阶可微，则可以使用牛顿法来求解函数的极值。牛顿法6.1一维搜索算法▌例6.4利用牛顿法确定函数f

(x)=0.5x2−sinx的极值点，初始值为x(0)=0.5，精度为ε=10−5，|x(k+1)−x(k)|<ε时停止循环。使用Julia语言在Syslab中编程，代码如下：clear()usingTySymbolicMath@variablesxf=0.5*x^2-sin(x)

#定义目标函数df=derivative(f,x)#定义目标函数的一阶导数d2f=derivative(df,x)#定义目标函数的二阶导数k0=0.5 #给定初始值epsilon=10^(-5)

#定义收敛要求dx=value(substitute(df,Dict(x=>k0)))d2x=value(substitute(d2f,Dict(x=>k0)))k1=k0-dx/d2xwhileabs(k0-k1)>epsilond2x=value(substitute(d2f,Dict(x=>k0)))globalf,df,d2f,k0,k1,epsilon,dx,d2xk0=k1dx=value(substitute(df,Dict(x=>k0)))k1=k0-dx/d2xendx_min=round(k1,digits=4)fmin=value(substitute(f,Dict(x=>x_min)))#获取函数的极小值f_min=round(fmin,digits=4)print("x_min=$x_min,f_min=$f_min")输出：x_min=0.7391,f_min=-0.40056.2线性优化算法线性规划研究的是一类在线性约束条件下求解线性目标函数极值的问题，即确定一组决策变量，使得目标函数取得极大值或极小值。线性规划是一类特殊的有约束优化问题，目标函数的极值通常是求取极小值。任何满足约束条件的点均称为可行点。在线性规划问题中，目标函数是线性的，可行点的集合由一组线性等式和/或不等式确定。6.2线性优化算法单纯形法矩阵：A是m×n实数矩阵，m<n，rand(A)=m。约束：不失一般性，假设b≥0。如果列向量b中的第i个元素是负数，那么在第i个约束方程两端同时乘以−1，就可以使方程的右端项大于零。转换为标准型：如果不等式约束条件为Ax≤b，那么也可以引入松弛向量，转换为标准型：Ax+Imy=b。标准型不等式形式引入剩余变量6.2线性优化算法令A的前m列是基向量，组成m×m非奇异矩阵B。A的非基列向量组成m×(n−m)矩阵D，价值系数可写为cT=[cBT

cDT]。标准型可表示为:xD=0如果xD=0，那么x=[xBT,xDT]T=[xBT,0]T=[B−1b,0]T是关于基B的基本可行解，相应的目标函数值为z0=cBTB−1b。xD≠0如果xD≠0，那么x=[xBT,xDT]T不是基本解。此时，x=[xBT,xDT]T可表示为xB=B−1b−B−1DxD，相应的目标函数值为z=cBTxB+cDTxD=cBTB−1b+(cDT−cBTB−1D)xD=z0+rDTxD向量rD=(cDT−cBTB−1D)T中的元素是非基变量的检验数。6.2线性优化算法构造单纯形法的增广矩阵形式：对增广矩阵进行初等行变换：在变换后的增广矩阵的最后一行中，cDT−cBTB−1D是检验数，−cBTB−1b是当前基本可行解下目标函数取负后的结果。6.2线性优化算法▌例6.5利用单纯形法求解下列函数的极值：为了利用单纯形法求解上述问题，先要引入松弛变量，把问题化成标准型：6.2线性优化算法编程运算所需的系数矩阵如下：更新增广矩阵标准型：6.2线性优化算法更新增广矩阵标准型：单纯形表的最后一行矩阵向量c全部非负，因此其对应的基本可行解是最优解，标准型目标函数的最小值为−136/3，原问题目标函数的最大值为136/3，最优基本可行解x=[4,0,8/3,0]。6.2线性优化算法使用Julia语言在Syslab中编程，代码如下：clear()#SimplexmethodE=[11011004.0;-14210106;1-13200112;-106-2-40000]#由于计算可能出现浮点数，因此将初始E矩阵设定为浮点矩阵flag=trueS=0#迭代次数Z=0.0m,n=size(E)BIndex=collect(n-m+1:n-1)#初始基变量索引（向量）ifn<merror("错误：变量个数不能少于约束数量（E的列数应大于或等于行数）")else

X=zeros(n-1)whileflagglobalE,m,n,BIndex,X,flag,Z,P,S,k1,k2ifminimum(E[m,:])>=0X.=0.0foriin1:m-1X[BIndex[i]]=E[i,n]endZ=dot(E[m,1:n-1],X)println("迭代次数为:$S")println("已找到最优解。")flag=falsebreakelse6.2线性优化算法S+=1_,k1=findmin(E[m,1:n-1])#确定进基向量位置

P=fill(Inf,m-1)foriin1:m-1ifE[i,k1]>0P[i]=E[i,n]/E[i,k1]endendifall(isinf.(P))error("问题无界：所有E[:,k1]<=0，无法继续迭代。")end输出：迭代次数为:2已找到最优解。最优解X=[4.0,0.0,2.667,0.0,0.0,4.667,0.0]最小目标值Z=-45.333_,k2=findmin(P)BIndex[k2]=k1E[k2,:]./=E[k2,k1]foriin1:mifi!=k2E[i,:]-=E[i,k1]*E[k2,:]endendendendendZ=-E[m,n]#目标函数值取负println("最优解X=",round.(X,digits=3))println("最小目标值Z=",round(Z,digits=3))6.2线性优化算法区别：单纯形法需要从一个标准的单纯形表开始运算，如果通过穷举法寻找初始基本可行解，则需要任选m个基列向量，并将单纯形表转换为标准型。整个穷举法过程可能需要尝试

次，比较烦琐。尤其在需要引入剩余变量时，没有明显的基本可行解，单纯形法无法开展。两阶段单纯形法通过引入人工变量求解一般形式的线性规划问题。基本思路第一阶段，引入人工变量，构造人工目标函数，找到基本可行解；第二阶段，将第一阶段得到的基本可行解作为初始条件，采用单纯形法求解原规划问题。两阶段单纯形法6.2线性优化算法▌例6.6求解下列函数的极值：首先引入剩余变量x3与x4，把问题化成标准型：该标准型没有明显的基本可行解，单纯形法无法开展，因此采用两阶段单纯形法。6.2线性优化算法第一阶段，引入人工变量x5,x6≥0，构造人工目标函数x5+x6，获得人工问题的单纯形表，经更新单纯形表的最后一行，可获得标准型如下：更新单纯形表：6.2线性优化算法对最后一行进行初等变换，使得与基向量对应的最后一行元素为0，可得：所有检验数都非负，最优解为x=[18/7,6/7,0,0]T，目标函数的最优值是54/7。6.2线性优化算法使用Julia语言在Syslab中编程，将单纯形法作为子函数，被两阶段单纯形法主函数调用。子函数Simplex_function.jl()的主体程序不变，主要改动如下如下：functionSimplex_function(E::Matrix{Float64})#改写为子函数格式#子函数主体部分省略returnX,Z,BIndexend输出：最优解X=[2.571,0.857,0.0,0.0,0.0,0.0]最小目标值Z=7.7146.2线性优化算法两阶段单纯形法的主函数代码如下：clear()#Two-stagesimplexmethodinclude("Simplex_function.jl")#调用自定义子函数globalE,m,n,BIndex#第一阶段#引入人工变量的增广矩阵EE=[42-101012.0;140-1016;0000110]m,n=size(E)foriin1:n-1localkifE[m,i]!=0forjin1:m-1ifE[j,i]!=0k=E[m,i]/E[j,i]E[m,:]=E[m,:]-k*E[j,:]BreakendendendendX,Z,BIndex=simplex_function(E)#第二阶段，更新E矩阵，去除引入的人工变量，更新目标函数A1=hcat(E[:,1:4],E[:,7])E=A1E[m,1]=2E[m,2]=3m,n=size(E)foriin1:m-1localkk=E[m,BIndex[i]]/E[BIndex[i],BIndex[i]]E[m,:]-=E[BIndex[i],:]*kendprintln("最优解X=",round.(X,digits=3))println("最小目标值Z=",round(-E[m,n],digits=3))6.3蚁群优化算法群体智能是一种模拟自然界群体行为的计算方法，借鉴了群体动物或昆虫在协作中展现出来的智能。在群体智能中，个体之间相互通信、相互协作，通过简单的规则和局部信息交流来实现整体上的智能行为。而优化算法则是一类用于解决最优化问题的数学方法，能够在大量搜索空间中找到最优解。6.3蚁群优化算法基本思路ACO是一种用来寻找优化路径的概率型算法。ACO源于对蚂蚁寻找食物路径行为的模拟。ACO通过模拟蚁群在环境中的寻找和选择过程来寻找最优解：蚂蚁在搜索过程中会释放信息素，其他蚂蚁根据信息素的浓度选择路径，最终形成一条最佳路径。ACO原理6.3蚁群优化算法步骤：初始化：首先确定蚂蚁的数量、问题的规模、信息素的初始浓度等参数，并将蚂蚁随机放置在各个起点上。蚂蚁移动：每只蚂蚁根据当前位置和信息素的浓度选择下一个要移动到的位置。蚂蚁选择路径的概率与路径上的信息素的浓度及启发式信息有关。启发式信息通常是根据问题的特点来定义的，如在旅行商问题中可以定义为两个城市相隔距离的倒数。信息素更新：当所有蚂蚁都完成一次遍历后，对路径上的信息素进行更新。信息素更新分为局部更新和全局更新。局部更新是指蚂蚁在走过的路径上实时更新信息素，使其浓度降低，以鼓励探索新的路径。全局更新是指在所有蚂蚁完成一次遍历后，根据蚂蚁找到的最优路径更新信息素，使最优路径上的信息素浓度提升。终止条件判断：判断是否达到预设的终止条件，如迭代次数达到上限、最优解不再变化等；如果未达到终止条件，则返回步骤（2）继续选代。6.3蚁群优化算法使用Julia语言在Syslab中编写高斯变异函数GaussMutation()，代码如下：（在局部更新中，使用高斯变异操作。高斯变异通过用正态分布随机数替换原有基因值来实现个体向量元素的随机扰动，其操作过程与均匀变异的操作过程类似。）#高斯变异函数functionGaussMutation(x,LB,UB,k,K)M=length(x)mx=similar(x)foriin1:M#变异强度随迭代衰减（前期探索，后期求精）

decay=1.0-(k-1)/K#从1.0衰减到0.0sigma=(UB[i]-LB[i])*0.15*decay#初始变异更强

mx[i]=x[i]+sigma*randn()mx[i]=clamp(mx[i],LB[i],UB[i])endreturnmxendlength(x)：返回元素数量；similar(x)：生成一个与x具有相同类型（这里指的是密集、稀疏等）的未初始化数组，且具有相同的元素类型和维数。randn()：生成一个随机数组，元素为标准正态分布，服从独立同分布。clamp(mx[i],LB[i],UB[i])：将元素mx的值限定在[LB,UB]内。6.3蚁群优化算法基本思路ACO参数初始化，主要包括迭代次数、蚁群规模、信息素蒸发系数、信息素增加强度、蚂蚁爬行速度、适应值等。设定被控对象及仿真参数，主要包括传递函数分子系数/分母系数、仿真步长等。计算当前代所有蚂蚁的位置、适应值，搜索获得最优蚂蚁位置及适应值，更新信息素和适应值。针对获得的当前代蚂蚁的最优位置和适应值，进行局部搜索，利用高斯变异生成新参数，计算变异后参数的适应值，并根据适应值更新蚂蚁位置。朝信息素的浓度最高的地方进行全局搜索，找到选中蚂蚁中信息素的浓度最高的位置，比较原位置和新位置的适应值，确定是否朝新位置移动。更新信息素，存储更新后的所有蚂蚁位置、适应值及最优位置和适应值。到达迭代停止条件，算法停止；否则回到步骤3，循环执行。基于ACO的PID参数整定6.3蚁群优化算法▌例6.7二阶系统的传递函数为：采用蚁ACO整定获得该系统的PID参数。使用Julia语言编写，设计主函数ACO_PID()、子函数PIDOBJ()和ACO_function()。子函数PIDOBJ()为评价函数，实现PID参数的适应值计算。针对二阶传递函数模型，以阶跃指令信号为例，评价函数主要包括响应时间和超调量等指标。ACO_function()实现ACO功能，需要调用6.3.1节的高斯变异函数GaussMutation()。整体算法需要调用多个Syslab的库文件，统一在主函数中进行引用。6.3蚁群优化算法使用Julia语言编写的主函数ACO_PID()的代码如下：#MainfunctionACO_PID()usingControlSystems #引用库函数usingTyPlotusingRandomusingTySignalProcessingusingXLSXinclude("ACO_function.jl") #调用自定义子函数include("GaussMutation.jl")include("PIDOBJ.jl")#参数设置

K=150 #迭代次数

N=50 #蚁群规模

Rho=0.8 #信息素蒸发系数

Q=1 #信息素增加强度Lambda=0.3 #蚂蚁爬行速度

#PID参数的上下界(Kp,Ki,Kd)LB=[0.1;0.0;0.0]#下界

UB=[100.0;20.0;20.0] #上界

#传递函数

Num=[133.0] #分子系数

Den=[1.0,25.0,0.0] #分母系数

Delay=0.0 #时间延迟

#仿真参数

ts=0.005 #仿真时间步长

StepNum=2000 #仿真总步数

SigType=1

#指令类型：1:阶跃,2:方波;3.正弦波6.3蚁群优化算法#PID输出限幅

PIDLB=-15.0 #输出下限

PIDUB=15.0 #输出上限

#调用ACO优化PID参数

BESTX,BESTY,ALLX,ALLY,k=ACO_function(K,N,Rho,Q,Lambda,LB,UB,Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)#获取最优PID参数

best_params=BESTX[k-1]Kp=best_params[1]Ki=best_params[2]Kd=best_params[3]#打印最优PID参数

println("最优PID参数:")println("Kp=",Kp)println("Ki=",Ki)println("Kd=",Kd)#验证最优参数J,u,yout,error,ITAE,OS,u_energy=PIDOBJ(Kp,Ki,Kd,Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)#打印评价指标

println("评价指标:")println("ITAE(时间加权绝对误差积分):",ITAE)println("超调量(OS):",round(OS*100,digits=2),"%")#转换为百分比

println("控制能量:",u_energy)println("综合评价函数值(J):",J)#绘制响应曲线

plt=TyPlot.plot(yout,label="系统输出")TyPlot.title("最优PID参数的系统响应")TyPlot.xlabel("时间步")TyPlot.ylabel("输出")display(plt)TyPlot.ylabel("输出")display(plt)6.3蚁群优化算法函数运行时，需要使用ControlSystems函数库文件。第一次使用时，需要首先在Julia终端输入“]”，进入pkg；然后输入addControlSystems即可安装。6.3蚁群优化算法PID参数的评价子函数PIDOBJ()的代码如下：#PID控制器的评价函数functionPIDOBJ(Kp,Ki,Kd,Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)#输入参数:#Kp,Ki,Kd-PID参数

#Num,Den-被控对象传递函数系数

#Delay-延迟时间

#输出参数:#J-评价函数值(越小越好)#u-控制输出

#yout-系统输出

#error-误差序列

#初始化

u=zeros(StepNum) #控制量

yout=zeros(StepNum)#系统输出

error=zeros(StepNum) #误差

integral=0.0 #积分项

error_prev=0.0 #上一时刻误差t=0:ts:(StepNum-1)*tsifSigType==1 #阶跃信号

r=ones(StepNum)elseifSigType==2 #方波信号

r=square(t.*2π*0.5).+1

else #正弦波信号

r=sin.(2π*0.5.*t)endsys=ControlSystems.tf(Num,Den)#被控对象模型（离散化）

ifDelay>0sys=sys*ControlSystems.delay(Delay)dsys=ControlSystems.c2d(sys,ts,:zoh)elsedsys=ControlSystems.c2d(sys,ts,:zoh)#零阶保持器离散化

end

6.3蚁群优化算法numd=num(dsys) #获取离散分子系数

dend=den(dsys) #获取离散分母系数

numd=numd[1] #转换为向量

n=length(numd) #分子系数长度

#仿真PID控制过程

forkin1:StepNum#当前误差

error[k]=r[k]-yout[k]#PID计算

proportional=Kp*error[k]integral+=Ki*error[k]*ts#积分项（累积误差×步长）

derivative=Kd*(error[k]-error_prev)/ts#微分项（误差变化率）

#控制量计算与限幅

u[k]=proportional+integral+derivativeu[k]=clamp(u[k],PIDLB,PIDUB)#输出限幅

#计算系统输出y_kifk<=n#当k小于分子系数的长度时，取u的前k个元素

u_slice=u[1:k]numd_slice=numd[1:k]#取numd的前k个元素

y_k=sum(numd_slice.*u_slice)/dend[1]else#当k大于或等于分子系数的长度时，取u的最近n个元素

u_slice=u[k-n+1:k]y_k=sum(numd.*u_slice)/dend[1]end#减去分母项的影响（离散系统递归计算）

foriin2:length(dend)ifk>=i-1y_k-=dend[i]/dend[1]*yout[k-i+2]endend#更新输出6.3蚁群优化算法ifk<StepNumyout[k+1]=y_kenderror_prev=error[k]#保存当前误差，用于下一时刻的微分计算

end#计算评价函数J(综合性能指标)t_vector=collect(t)ITAE=sum(t_vector.*abs.(error))#时间加权绝对误差积分

#计算超调量

y_max,idx=findmax(yout)OS=ifSigType==1&&r[1]!=0#阶跃响应的超调量计算

max(0.0,(y_max-r[1])/r[1])#超调量（百分比）

else0.0end#控制能量u_energy=sum(u.^2)*ts#综合评价函数(权重可调整)ifSigTypein[2,3]#方波/正弦波

#计算误差变化率

error_diff=diff(error)#误差的一阶差分

error_diff=vcat(error_diff,0.0)#保持长度一致

error_diff_penalty=sum(abs.(error_diff))*ts#误差变化率积分

#增大ITAE权重，同时惩罚误差变化

J=(10*ITAE+2*error_diff_penalty+0.05*u_energy)/10else#阶跃

J=ITAE+100*OS+0.1*u_energyEndreturnJ,u,yout,error,ITAE,OS,u_energyend6.3蚁群优化算法ACO子函数ACO_function()的代码如下：functionACO_function(K,N,Rho,Q,Lambda,LB,UB,Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)#输出参数如下

#Best_info：每一代的最优蚂蚁参数信息

#Best_fit：每一代的最优评价函数值

#All_pos：每一代所有蚂蚁的位置

#All_fit：每一代所有蚂蚁的评价函数值

M=length(LB)

#决策变量的个数（此处为3，对应Kp、Ki、Kd）

#蚁群位置初始化（每列代表一只蚂蚁的PID参数）

X=zeros(M,N)foriin1:Mx=LB[i].+(UB[i].-LB[i]).*rand(1,N)#在[LB,UB]范围内随机生成参数

X[i,:]=xEnd#输出变量初始化

All_pos=Vector{Matrix{Float64}}(undef,K)#存储每一代所有蚂蚁的位置（M×N矩阵）

All_fit=zeros(K,N)

#存储每一代所有蚂蚁的评价函数值

Best_info=Vector{Vector{Float64}}(undef,K)

#存储每一代的最优蚂蚁参数（M×1向量）

Best_fit=zeros(K)

#存储每一代的最优评价函数值

k=1#迭代计数器初始化

Tau=ones(N) #信息素初始化（每只蚂蚁对应一个信息素值）

Fit=zeros(N) #适应值初始化（存储当前代蚂蚁的评价函数值）6.3蚁群优化算法stop_size=10 #迭代结果相同的次数，迭代停止条件之一

recent_best_J=Float64[] #存储最近的最优J值（动态更新，最多保留stop_size个）

stop_flag=false

#停止迭代的标志（true=满足条件，跳出循环）

#迭代过程

whilek<=K&&!stop_flag#计算当前代所有蚂蚁的初始适应值

All_fit=zeros(N)fornin1:Nx=X[:,n]#第n只蚂蚁的参数（Kp,Ki,Kd）

#调用PIDOBJ()函数计算评价函数值JJ,u,yout,error=PIDOBJ(x[1],x[2],x[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)All_fit[n]=Jend

#确定当前代的最优蚂蚁位置（POS）min_fit=minimum(All_fit)#寻找最小的评价函数值（J越小越优）

temp_pos=findall(All_fit.==min_fit)POS=temp_pos[1]

#最优蚂蚁的索引（若有多个最优解，则取第一个）

#蚂蚁随机探路（局部搜索：高斯变异）

fornin1:Nifn!=POS#跳过当前最优蚂蚁

x=X[:,n]#当前蚂蚁的参数

#计算原参数的适应值

J,u,yout,error=PIDOBJ(x[1],x[2],x[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fx=J#高斯变异生成新参数

mx=GaussMutation(x,LB,UB,k,K)

6.3蚁群优化算法#计算变异后参数的适应值

J,u,yout,error=PIDOBJ(mx[1],mx[2],mx[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fmx=J#根据适应值更新蚂蚁位置

ifFmx<Fx#变异后更优，直接更新

X[:,n]=mxFit[n]=Fmxelseifrand()>1-(1/sqrt(k))#一定概率接受新解（概率随迭代降低）

Fit[n]=Fmxelse#保留原参数

X[:,n]=xFit[n]=Fxendendendfornin1:N#二次局部搜索（增强局部探索能力）

ifn!=POS#跳过当前最优蚂蚁

x=X[:,n]#当前蚂蚁的参数

#计算原参数的适应值

J,u,yout,error=PIDOBJ(x[1],x[2],x[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fx=Jmx=GaussMutation(x,LB,UB,k,K)

#再次进行高斯变异

#计算变异后参数的适应值

J,u,yout,error=PIDOBJ(mx[1],mx[2],mx[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fmx=JifFmx<Fx#变异后更优，更新适应值#更新规则（与第一次局部搜索略有差异）

Fit[n]=Fmx6.3蚁群优化算法elseifrand()>1-(1/sqrt(k))#概率接受新解并更新位置

X[:,n]=mxFit[n]=Fmxelse#保留原参数

X[:,n]=xFit[n]=Fxendendendfornin1:N#朝信息素的浓度最高的地方移动（全局搜索）

ifn!=POS#跳过当前最优蚂蚁

x=X[:,n]#当前蚂蚁的参数

r=(K+k)/(K+K)#选择比例随迭代增加#选择部分蚂蚁的信息素作为参考

p=shuffle(1:N)#随机排列蚂蚁索引

t=ceil(Int,r*N)#选择t只蚂蚁

pos=p[1:t]TempTau=Tau[pos]#选中蚂蚁的信息素

maxTempTau=maximum(TempTau)#找到选中蚂蚁中信息素的浓度最高的位置

pos2=findall(TempTau.==maxTempTau)pos3=pos[pos2[1]]x2=X[:,pos3[1]]#高信息素浓度蚂蚁的参数

x3=(1-Lambda)*x+Lambda*x2

#向高信息素浓度位置移动（加权平均）

#计算原位置和新位置的适应值

J,u,yout,error=PIDOBJ(x[1],x[2],x[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fx=J

J,u,yout,error=PIDOBJ(x3[1],x3[2],x3[3],Num,Den,Delay,ts,StepNum,SigType,PIDLB,PIDUB)Fx3=JifFx3<Fx#移动后更优，更新位置

#根据适应值更新位置6.3蚁群优化算法X[:,n]=x3Fit[n]=Fx3elseifrand()>1-(1/sqrt(k))#概率接受新位置

X[:,n]=x3Fit[n]=Fx3else#保留原位置

X[:,n]=xFit[n]=FxendendendTau=Tau.*(1-Rho)#信息素蒸发#更新信息素

max_fit=maximum(Fit)min_fit=minimum(Fit)ifmax_fit≈min_fit#优质解（Fit越小）获得更多信息素奖励

DeltaTau=zeros(N)elseDeltaTau=(max_fit.-Fit)./(max_fit-min_fit)#归一化处理

endTau=Tau.+Q.*DeltaTau#信息素更新

#记录当前代结果

All_pos[k]=copy(X) #存储当前代所有蚂蚁位置

All_fit[k,:]=Fit

#存储当前代所有适应值

currcent_min_fit=minimum(Fit) #当前代最优适应值

pos4=findall(Fit.==currcent_min_fit)Best_info[k]=X[:,pos4[1]] #当前代最优参数信息

Best_val[k]=currcent_min_fit #当前代最优适应值

push!(recent_best_J,currcent_min_fit)6.3蚁群优化算法iflength(recent_best_J)>stop_sizepopfirst!(recent_best_J)endiflength(recent_best_J)==stop_size#检查所有元素是否与第一个元素相等（浮点比较用isapprox）

all_same=all(isapprox.(recent_best_J,recent_best_J[1]))ifall_samestop_flag=true#满足条件，设置停止标志

println("最近",stop_size,"次最优J值全部相同")println("迭代提前终止，最终迭代次数：",k)println("最终最优J值：",round(currcent_min_fit,digits=4))endendif(k%10==0)println("迭代次数：",k,"，当前最优J：",round(currcent_min_fit,digits=4))endk+=1#迭代计数器加1endresults=Matrix{Any}(undef,k,5)fill!(results,missing)#初始化所有元素为missing，避免未定义引用

results[1,1]="迭代次数"results[1,2]="Kp"results[1,3]="Ki"results[1,4]="Kd"results[1,5]="最优评价函数值"ifSys.iswindows()desktop_path=joinpath(homedir(),"Desktop/ACO")end6.3蚁群优化算法excel_filename="ACO_Best_Parameters.xlsx"full_excel_path=joinpath(desktop_path,excel_filename)max_iter=min(k-1,K)foriin1:max_iterresults[i+1,1]=i#迭代次数

results[i+1,2]=Best_info[i][1]#Kp参数

results[i+1,3]=Best_info[i][2]#Ki参数

results[i+1,4]=Best_info[i][3]#Kd参数

results[i+1,5]=Best_val[i]#最优评价函数值

endXLSX.openxlsx(full_excel_path,mode="w")doxf#创建Excel文件并写入数据

sheet=xf[1]#获取默认工作表

sheet["A1"]=results #从A1单元格开始写入所有数据

endprintln("文件路径:",full_excel_path)returnBest_info,Best_val,All_pos,All_val,kend6.3蚁群优化算法输出：julia>正在运行ACO_PID.jl迭代次数：10，当前最优J：0.7251迭代次数：20，当前最优J：0.6926迭代次数：30，当前最优J：0.5983迭代次数：40，当前最优J：0.5526迭代次数：50，当前最优J：0.449最近10次最优J值全部相同迭代提前终止，最终迭代次数：59最终最优J值：0.449在该子函数中，通过引用库函数XLSX，将迭代产生的PID参数历史数据记录在Excel文件中。信号类型选择1：SigType=1，运行主函数ACO_PID.jl()。