CN110162841B - 一种引入三维稳定性约束的铣削加工多目标优化决策方法 - Google Patents

一种引入三维稳定性约束的铣削加工多目标优化决策方法 Download PDF

Info

Publication number
CN110162841B
CN110162841B CN201910342927.6A CN201910342927A CN110162841B CN 110162841 B CN110162841 B CN 110162841B CN 201910342927 A CN201910342927 A CN 201910342927A CN 110162841 B CN110162841 B CN 110162841B
Authority
CN
China
Prior art keywords
firefly
machine tool
milling
constraint
fireflies
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201910342927.6A
Other languages
English (en)
Other versions
CN110162841A (zh
Inventor
杨文安
蔡旭林
黄超
郭栋琨
杨博
郭宇
廖文和
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University of Aeronautics and Astronautics
Original Assignee
Nanjing University of Aeronautics and Astronautics
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University of Aeronautics and Astronautics filed Critical Nanjing University of Aeronautics and Astronautics
Priority to CN201910342927.6A priority Critical patent/CN110162841B/zh
Publication of CN110162841A publication Critical patent/CN110162841A/zh
Application granted granted Critical
Publication of CN110162841B publication Critical patent/CN110162841B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Biophysics (AREA)
  • Artificial Intelligence (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Computational Linguistics (AREA)
  • Biomedical Technology (AREA)
  • Mathematical Physics (AREA)
  • Software Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Numerical Control (AREA)

Abstract

本发明公开了一种引入三维稳定性约束的铣削加工多目标优化决策方法，包括：构建关于铣削主轴转速、轴向切深以及径向切深的三维稳定性模型；建立以高效率、低成本和高利润为目标函数，以机床性能、刀具性能、工件性能、稳定性要求以及加工要求为约束条件的铣削参数优化模型；利用基于聚类原理的多目标萤火虫优化算法求解铣削参数优化模型，得到多组最优工艺参数。本发明针对铣削加工工艺参数的多目标优化提出了更为合理的优化决策方法，为多次走刀高速铣削工艺参数的优化选择提供了一种有效且实用的技术。

Description

一种引入三维稳定性约束的铣削加工多目标优化决策方法
技术领域
本发明属于机械加工工艺领域，尤其涉及一种引入三维稳定性约束的铣削加工多目标优化决策方法。
背景技术
数控技术作为现代先进制造装备行业中的核心技术，其切削参数的合理选择极大的影响了加工效率、加工质量以及加工成本，因此国内外的学者针对铣削加工工艺参数优化进行了大量研究。
国内，公开号为CN201410487219.9的发明专利公开了一种机械加工能量效率的切削参数优化方法，该方法未对径向切深进行优化，且未考虑铣削颤振稳定性约束的影响，因此用所获得的工艺参数进行加工极有可能导致振纹及崩刃等问题。公开号为CN201310109671.7的发明专利公开了一种基于特征的飞机结构件切削参数优化方法，该方法不能保证铣削的稳定性，且该方法在实际生产中操作过于繁琐，鲁棒性无法预估，因而工程应用意义不大；除此之外，该专利中对于优化的多个不同量纲的目标采用加权相加的方法转化为单目标进行处理，然而在相加之前需要对目标函数进行标定，但目标值的最大、最小边界不可预测，因此该方法是不合理的，且这样不能同时得到多组具有不同特性的非劣最优解，无法真正发挥多目标优化的优势。公开号为CN201410201398.5的发明专利公开了一种铣削工艺参数优化方法，该方法未对径向切深进行优化，且该专利同样采用不同量纲的目标函数加权相加的方法转化为单目标进行处理。胡瑞飞等(胡瑞飞,殷国富,殷鸣.切削稳定性约束下的铣削参数优化技术研究[J].机械工程学报,2017,53(5):190-198.)提出的该方法未考虑加工粗糙度，主轴功率等约束；所用的铣削稳定性约束模型虽然同时考虑了主轴转速，轴向切深以及径向切深，但是这里的约束是多个二维稳定性约束所构成的伪三维稳定性约束，并不能找到真正意义上的最优解。
国外，Pawar P J等(Pawar P J,Rao R V.Parameter optimization ofmachining processes using teaching–learning-based optimization algorithm[J].The International Journal of Advanced Manufacturing Technology,2013,67(5-8):995-1006.)提出的方法未对径向切深进行优化，且未考虑加工粗糙度及颤振稳定性的约束，所得参数实用性不强。Yildiz A R(Yildiz A R.A new hybrid differentialevolution algorithm for the selection of optimal machining parameters inmilling operations[J].Applied Soft Computing,2013,13(3):1561-1566.)提出的方法未对径向切深进行优化，且未考虑颤振稳定性的约束，因此所得参数实用性不强。RashidM.F.F Ab等(Ab Rashid M F F,Harun W S W,Ghani S A C,et al.Optimization ofMulti-Pass Pocket Milling Parameter Using Ant Colony Optimization[J].AdvancedMaterials Research,2014,1043：65-70.)提出的方法未考虑颤振稳定性、刀具寿命、机床功率等约束；并且分配轴向切深的方法为粗略的平均分，这样处理在实际生产中并不合理。KhalilpourazariS等(Khalilpourazari S,Khalilpourazary S.SCWOA:An efficienthybrid algorithm for parameter optimization of multi-pass milling process[J]，Journal of Industrial and Production Engineering，2018,35(3)135-147)提出的方法未考虑颤振稳定性约束；分配轴向切深的方法为粗略的平均分未考虑粗加工及精加工，因此无法有效的保证加工的稳定性及产品的表面粗糙度。
关于铣削三维稳定性预测方法，Vincent Thevenot等(Vincent Thevenot,LionelArnaud,Gilles Dessein,Gilles Cazenave-Larroche.Integration of dynamicbehaviour variations in the stability lobes method.Int.J.Adv.Manuf.Technol,2006(27):638-644)提出的该方法的局限性体现在以下几个方面：他所建立的分析模型采用直角切削模型建立了特征方程，与实际切削情况有一定的差别，且未考虑材料去除对于颤振稳定性的影响，更重要的是，没有考虑径向切深和稳定性的关系。U.Bravo等(UBravo,OAltuzarra,LN López de Lacalle.Stability limits of milling considering theflexibility of the workpiece and the machine.International Journal of MachineTools&Manufacture,2005(45):1669-1680)提出的方法中只是考虑了有限的加工阶段，且只是通过有限元模拟出不同加工阶段的三维稳定性图，局限性较大。Tony L.Schmitz等基于ZOA法(零阶法，zero order approach)构建了三维稳定性图，但是ZOA法在小径向切深时预测不准的问题也同样导致了该方法不适用于精加工等去除量较小的加工阶段。
现有技术中萤火虫算法只有个别用于多目标问题的优化中，但最终所获得的解的多样性及延展性均不够理想，且未见将萤火虫算法应用于切削参数优化问题中的研究；目前针对多目标优化的研究均仅为获得大量Pareto前沿解，而并未考虑大量的输出并不利于工艺师进行决策。
发明内容
发明目的：为解决上述现有技术中约束条件考虑的不够完善、局限性大、预测不准准确、操作繁琐等问题，本发明提供一种引入三维稳定性约束的铣削加工多目标优化决策方法。
技术方案：本发明提供一种引入三维稳定性约束的铣削加工多目标优化决策方法，包括如下步骤：
步骤1：基于二阶全离散法，在主轴转速Ω、轴向切深d的铣削稳定性预测模型中引入径向切深a，构建关于主轴转速Ω、轴向切深d以及径向切深a的铣削三维稳定性预测模型；
步骤2：建立多目标多约束的多次走刀铣削参数优化模型，所述该模型包括以高效率、低成本以及高利润为目标的多目标函数；和根据数控机床性能、刀具性能、工件性能、加工要求以及步骤1中的三维稳定性模型建立的多约束条件；
步骤3：利用基于聚类原理的多目标萤火虫优化算法求解多目标多约束的多次走刀铣削参数优化模型，得到多组具有代表性的切削参数，即多组具有代表性的Pareto前沿解，所述基于聚类原理的多目标萤火虫优化算法包括：基于云模型的坐标摄动机制、基于拥挤度的自适应吸引力机制、基于拥挤度的最优萤火虫引导机制、基于聚类原理的外部存档更新消减机制；
步骤4：利用伪权重估计法，计算每个具有代表性的前沿解其不同目标函数的权重值。
进一步的，所述径向切深a的表达式为：
Figure BDA0002041337110000031
其中s为无穷小量，ε为径向切深被平均划分为ε段，D为铣刀直径。
进一步的，所述步骤2中的多目标函数包括：铣削耗时的目标函数、铣削成本的目标函数、铣削利润的目标函数；
所述铣削耗时的目标函数如下所示：
Figure BDA0002041337110000032
其中，D为铣刀直径；Ωri为第i次粗加工时机床的主轴转速，Ωs为精加工时机床的主轴转速，n为粗加工的次数，z为铣刀齿数，fri为第i次粗加工时候每齿进给量；fs为精加工时每齿进给量；h1以及h2为与走刀长度以及切入、切出时间相关的常数；Trc为换刀时间；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命，Lri为第i次粗加工时的走刀长度，LS为精加工时的走刀长度；
Figure BDA0002041337110000041
Figure BDA0002041337110000042
Figure BDA0002041337110000043
Figure BDA0002041337110000044
其中dri为第i次粗加工时机床的轴向切深；ari为第i次粗加工时机床径向切深，Ωs为精加工时机床的主轴转速；as为精加工时机床的径向切；Cv,Kv,xv,yv,sv,qv,pv和l是与刀具和工件材料相关的常数；GInt()和SInt()分别表示取最大整数和最小整数，W和L分别为工件的宽度以及长度，e’是为避免撞刀而设置的任意距离；
所述铣削成本的目标函数如下所示：
Figure BDA0002041337110000045
其中：k0为单位时间的人工及管理成本，kt为单位时间的加权刀具成本，Tp为加工准备时间；
所述铣削利润的目标函数如下所示：
Figure BDA0002041337110000046
其中：Sp为产品的市场价格，Cmat为原材料成本。
进一步的，所述计算具有代表性的前沿解其不同目标函数的权重值的具体方法如下所示：
Figure BDA0002041337110000051
其中，wy为某个前沿解的第y个目标函数的权重值，fy,max为所有前沿解的第y个函数的最大值；fy为该某个前沿解第y个目标函数的值；fy,min为所有前沿解的第y个函数的最小值；Emax为所有目标函数的总个数，fE,max为所有前沿解的第E个函数的最大值，fE为该某个前沿解的第E个目标函数的值，fE,min为所有前沿解的第E个函数的最小值。
进一步的，所述多约束条件包括：主轴转速约束条件、进给速度约束条件、轴向切深约束条件、径向切深约束条件、切削力约束条件、切削力矩约束条件、主轴功率约束条件、表面粗糙度约束条件、刀具寿命约束条件和三维稳定性模型构成的约束条件；
所述主轴转速约束条件为：
Ωmin＜Ωri＜Ωmax；Ωmin＜Ωs＜Ωmax；
其中，Ωmax以及Ωmin为机床的最大、最小主轴转速；Ωri为第i次粗加工时机床的主轴转速；Ωs为精加工时机床的主轴转速；
所述进给速度约束条件为：
fmin＜fri＜fmax；fmin＜fs＜fmax；
其中，fmax以及fmin为机床的最大、最小每齿进给量即进给速度；fri为第i次粗加工时机床的每齿进给量；fs为精加工时机床的每齿进给量；
所述轴向切深约束条件为：
dmin＜dri＜dmax；dmin＜ds＜dmax；
Figure BDA0002041337110000052
其中，dmax以及dmin为机床的最大、最小轴向切深；dri为第i次粗加工时机床的轴向切深；ds为精加工时机床的轴向切深；dt为总铣削深度；n为粗加工的总次数；
所述径向切深约束条件为：
amin＜ari＜amax；amin＜as＜amax；
其中，amax以及amin为机床的最大、最小径向切深；ari为第i次粗加工时机床的径向切深；as为精加工时机床的径向切深；
所示切削力约束条件为：
Figure BDA0002041337110000061
其中，Cu,Ku,pu,qu,su,xu以及yu为常数，Fmin以及Fmax为机床的最大、最小主切削力；Fri为第i次粗加工时机床的主切削力；Fs为精加工时机床的主切削力；
所述切削力矩约束条件为：
Figure BDA0002041337110000062
其中：TMmax为机床的最大主轴转矩；TMri为第i次粗加工时机床的主轴转矩；TMs为精加工时机床的主轴转矩；
所述主轴功率约束条件为：
Figure BDA0002041337110000063
其中：pmax为机床的最大主轴功率；Pri为第i次粗加工时机床的主轴功率；Ps为精加工时机床的主轴功率；
所述表面粗糙度约束条件为：
Figure BDA0002041337110000064
其中，re为刀尖圆角半径，Rrmax以及Rsmax分别为粗加工和精加工最大表面粗糙度；Rri为第i次粗加工时表面粗糙度；Rs为精加工时表面粗糙度；
所述刀具寿命约束条件为：
tri≥TR，ts≥TR；
其中，TR为使用者所要求的刀具最短加工寿命；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命。
进一步的，所述得到若干个具有代表性的解的具体步骤如下所示：
步骤3.1：对每只萤火虫进行编码，并根据所设定编码方案和编码中的上边界L和下边界U初始化每只萤火虫的坐标，具体编码方式如下所示：
pA＝(n，dr，dh，ds，ar，ah，as，fr，fh，fs，Ωr，Ωh，Ωs)
L＝(nmin，drmin，dhmin，dsmin，armin，ahmin，asmin，frmin，fhmin，fsmin，Ωrmin，Ωhmin，Ωsmin)
U＝(nmax，drmax，dhmax，dsmax，armax，ahmax，asmax，frmax，fhmax，fsmax，Ωrmax，Ωhmax，Ωsmax)
其中，pA为第A只萤火虫的编码，n为粗加工的次数，dh为半精加工时的轴向切深，ah为半精加工时径向切深，fh为半径加工时的进给速度；Ωh为半精加工时主轴转速；
步骤3.2：计算初始种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度G(pA)；所述约束违背程度G(pA)的具体计算方法如下所示：
Figure BDA0002041337110000071
GC(pA)＝max{0，gC(pA)}
其中，G(pA)为第A个萤火虫的总约束违背程度，GC(pA)为第A个萤火虫的第C个约束条件的违背程度；CMAX为约束条件的总数量；若在计算时存在等式约束H(p)，则将其转变为不等式约束：|H(p)|-δ≤0,其中δ＝1.0E-004；
步骤3.3：根据选择Pareto前沿解的标准在当前萤火虫种群中获取Pareto前沿解，将所述Pareto前沿解放入外部存档中，如果外部存档超过预设的范围，则利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤3.4：在当前Pareto前沿解中选择最优萤火虫；具体选择方法为：将Pareto前沿解中的所有萤火虫按照拥挤度的升序排序，并在排序的顶部选择一定比例的萤火虫作为最优萤火虫候选，将候选的萤火虫按照每个目标函数值的大小进行排序，选择任意两个毗邻的萤火虫，以这两个萤火虫坐标连线的中点作为最优萤火虫的坐标；
步骤3.5：比较当前种群中任意两个萤火虫A与H的支配关系，若A支配H，则转步骤3.6；否则转步骤3.7；
步骤3.6：利用基于云模型的坐标摄动机制处理萤火虫A的坐标，利用基于拥挤度的自适应吸引力机制更新萤火虫H的坐标；并转向步骤3.8；
步骤3.7：利用基于云模型的坐标摄动机制处理最优萤火虫的坐标，利用基于拥挤度的最优萤火虫引导机制更新萤火虫H坐标；
步骤3.8：每个萤火虫的坐标均被移动过之后，按照一定的比例选取萤火虫，再次利用基于云模型的坐标摄动机制对选取的萤火虫的坐标进行处理；计算当前种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；并转步骤3.9和步骤3.10；
步骤3.9：根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，将其暂时存入外部存档，并转向步骤3.12；
步骤3.10：利用基于云模型的坐标摄动机制对步骤3.8获得的种群中的每只萤火虫的坐标进行处理，并计算处理后的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；
步骤3.11：再次根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，并将其暂时存入外部存档；
步骤3.12：对外部存档中的所有的Pareto前沿解按照选择Pareto前沿解的标准，重新获取Pareto前沿解，并剔除非前沿解；
步骤3.13：如果外部存档超过预设的范围，利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤3.14：IterNum＝IterNum+1；IterNum＞IterNummax？如果是转步骤3.15，否则转步骤3.4；
步骤3.15：输出外部存档中的所有萤火虫，即多组具有代表性的Pareto前沿解。
进一步的，在每次迭代计算之前基于聚类原理的多目标萤火虫优化算法会根据总切削深度dt，针对多次走刀高速铣削中的粗加工、半精加工以及精加工的轴向切削的深度进行划分，具体划分方案为：
Figure BDA0002041337110000081
则前n次走刀为粗加工,第n+1次走刀为精加工,如果n不是整数，则粗加工次数
为；
Figure BDA0002041337110000082
dh＝dt-ds-n*dr
则前n次走刀为粗加工，第n+1次走刀为半精加工；第n+2次走刀为精加工。
进一步的，所述选择Pareto前沿解的标准为：将每个萤火虫与其他萤火虫依次组成一对，如果两个萤火虫均为可行解即G(pA)＝0，在该对中选择占主导地位的萤火虫为该对的前沿解，如果两者均不占主导地位，则随机选择一个萤火虫作为该对的前沿解，如果一对中只有一个可行解则选择可行解为该对的前沿解，如果两者均不是可行解，则选择G(pA)最小的萤火虫作为该对的前沿解，所述占主导地位的萤火虫为Ttotal值小、Ctotal值小、ptotal值大的萤火虫。
进一步的，所述基于拥挤度的自适应吸引力机制为：
Figure BDA0002041337110000091
其中，
Figure BDA0002041337110000092
为本次迭代计算中萤火虫H被基于拥挤度的自适应吸引力机制处理后的坐标；
Figure BDA0002041337110000093
本次迭代计算中萤火虫H被处理前的坐标；
Figure BDA0002041337110000094
为萤火虫A的坐标；
Figure BDA0002041337110000095
为萤火A对萤火H的吸引度，β0为r＝0时萤火虫的吸引力；γ为光吸收系数；α′为随机化参数；
Figure BDA0002041337110000096
为与
Figure BDA0002041337110000097
大小成正相关的随机数向量；Rand(SizeVar)为[0,1]之间的SizeVar大小的随机数矩阵，SizeVar为设计变量的规模，Distance(·)表示求解的距离；CrowdH为萤火虫H的拥挤度；
所述基于拥挤度的最优萤火虫引导机制如下所示：
Figure BDA0002041337110000098
其中，
Figure BDA0002041337110000099
为第IterNum次迭代计算中萤火H被基于拥挤度的最优萤火虫引导机制处理后的坐标，
Figure BDA00020413371100000910
为萤火H被处理前的坐标，
Figure BDA00020413371100000911
为最优萤火虫对萤火虫H的吸引力，
Figure BDA00020413371100000912
为最优萤火虫的坐标；
所述基于云模型的坐标摄动机制的公式如下所示：
Figure BDA0002041337110000101
其中，En以及He分别表示随着迭代次数的增加而减小的熵和超熵，
Figure BDA0002041337110000102
表示被基于云模型的坐标摄动机制处理前萤火中A的坐标，
Figure BDA0002041337110000103
表示被处理后的萤虫A的坐标；
Figure BDA0002041337110000104
表示被基于云模型的坐标摄动机制处理前最优萤火O的坐标；
Figure BDA0002041337110000105
表示被处理后最优萤火O的坐标。
进一步的,所述利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减的具体操作如下所示：
步骤A：初始化集群集Q，该集群集中包括本次迭代计算得到的所有的Pareto前沿解；每个Pareto前沿解构成一个独立的集群；
步骤B：计算集群集中，每个集群与其他集群之间的距离d的计算公式如下所示，
Figure BDA0002041337110000106
其中：||·||代表为欧式距离；I1为集群Cδ中的萤火虫，I2为集群Cβ中的萤火虫，Cδ∈Q，Cβ∈Q；
步骤C：选择两个距离最小的集群，并用这两个集群的并集代替这两个集群，则集群集Q的变化如下所示：
Q＝Q-{Cδ，Cβ}∪{Cδ∪Cβ}
步骤D：并计算该集群集Q中集群的数量是否减少到预先指定的数量，如果是则转步骤E；否则转步骤B；
步骤E：最终得到指定数量的集群，判断每个集群中的萤火虫数量是否为1，如果集群中萤火虫的数量为2，则选择拥挤度大的萤火为这个集群的代表，如果数量超过2选择与该集群中其他萤火虫平均距离最小的萤火虫作为该集群的代表。
有益效果：本发明通过搭建三维稳定性模型、铣削参数优化模型以及针对多次走刀、多目标铣削参数优化的多目标萤火虫算法，实现了在机床性能、工件性能、刀具性能以及关于主轴转速、轴向切深和径向切深的三维颤振稳定性约束下的高速铣削每次走刀的主轴转速、进给量、轴向切深、径向切深的优化，以获得真正意义上的高效率，低成本、高利润率的具有代表性的多组最优解；且该多组最优解具有多样性、延展性以及较好的优化程度，除此之外，提出了一种针对多目标决策的伪权重计算方法对所获得的具有代表性的多组Pareto前沿解进行处理，以利于工艺师进行决策。
附图说明
图1为本发明的基于聚类原理的多目标萤火虫优化算法的流程图；
图2为本发明的多次走刀逆铣的走刀路线图；
图3为本发明的三维稳定性模型中的参数所绘制的三维稳定性叶瓣图；
图4为本发明的拥挤度计算方法；
图5为本发明的铣削参数优化模型的求解结果以及决策结果，其中(a)为500组Pareto前沿解；(b)为利用更新消减机制后具有代表性的20组前沿解；
图6为本发明的算法在测试函数中优化的结果；其中(a)是算法在测试函数CONSTR中的优化结果；(b)是算法在测试函数TNK中的优化的结果；(c)是算法在测试函数KITA中的优化的结果；(d)是算法在测试函数SRN中的优化的结果；(e)是算法在测试函数BNH中的优化的结果；(f)是算法在测试函数OSY中的优化的结果。
具体实施方式
构成本发明的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
本实施例提供一种引入三维稳定性约束的铣削加工多目标优化决策方法，包括如下步骤
步骤一：通过锤击激振法，对铣削系统的模态参数进行识别获取；
步骤二：并利用模态参数和二阶全离散法，构建关于主轴转速Ω、轴向切深d以及径向切深a的铣削三维稳定性模型；
步骤三：根据数控机床铣削过程建立以高效、低成本以及高利润为对象的铣削工艺参数优化模型的目标函数；
步骤四：根据数控机床性能、刀具性能、工件性能、加工要求以及步骤二中的三维稳定性模型建立铣削工艺参数优化模型的约束条件；
步骤五：根据实际生产需求，针对多次走刀高速铣削进行切削深度的划分，构建了一种基于聚类原理的用于求解大量约束条件下的多目标铣削工艺参数优化的多目标萤火虫优化算法，该算法包括基于云模型的坐标摄动机制、基于拥挤度的自适应吸引力机制、基于拥挤度的最优萤火虫引导机制、基于聚类原理的外部存档更新消减机制；通过该算法使得最终所获得的Pareto前沿解优化程度高、多样性高、延展性好以及代表性强；
步骤六：利用伪权重估计法，计算每个具有代表性的前沿解其不同目标函数的权重值。
步骤二所述利用模态参数和二阶全离散法，构建关于主轴转速Ω、轴向切深d以及径向切深a的铣削三维稳定性模型，具体包括如下步骤：
步骤(2-1)构建考虑再生颤振的铣削动力学模型：
Figure BDA0002041337110000121
其中ζ为相对阻尼，wn为固有频率，d为轴向切深，mt为模态质量，T为颤振周期，h(t)为切削力系数，x(t)为刀具的位移。
步骤(2-2)构建基于二阶全离散法的铣削三维稳定性模型：
考虑再生效应的铣削动力学模型用以下状态空间形式描述：
Figure BDA0002041337110000122
其中：
Figure BDA0002041337110000123
将周期T离散为m段,则步长
Figure BDA0002041337110000124
将
Figure BDA0002041337110000125
的表达式在步长τ的范围内进行积分，得到如下公式：
Figure BDA0002041337110000131
Figure BDA0002041337110000132
对x(t)、x(kτ+τ)进行线性插值以及拉格朗日变换得到：
Figure BDA0002041337110000133
x(kτ+τ-ζ-T)＝xk+1-m+ζ(xk-m-xk+1-m)/τ
Figure BDA0002041337110000134
并将A(kτ+τ-ξ)、x(kτ+τ-ξ-T)、x(kτ+τ-ξ)的表达式代入x(t)、x(kτ+τ)得：
xk+1＝Fk+1xk+1+(F0+F0,k)xk+Fk-1xk-1-Fk+1-mxk+1-m-Fk-mxk-m
其中：
Figure BDA0002041337110000135
F0＝Φ4
Figure BDA0002041337110000136
Figure BDA0002041337110000137
Figure BDA0002041337110000138
Figure BDA0002041337110000139
其中：
Figure BDA00020413371100001310
Figure BDA00020413371100001311
若矩阵[I-Fk+1]是非奇异的，则离散映射为:
Xk+1＝DkXk
其中：
Xk＝col(xk,xk-1,...,xk+1-m,xk-m)
Figure BDA0002041337110000141
利用离散映射序列Dk建立周期时间间隔上的转移矩阵：
Xm＝Dm-1Dm-2...D1D0X0
根据Floquet理论可以得到转移矩阵的所有特征值的模量，取模量等于1的特征值建立主轴转速Ω和轴向切深d的稳定性边界，即稳定性模型。
三维稳定性模型的第三个变量是径向切深a。为了得到理想的三维稳定性模型，定义无穷小量s，将径向切深的范围平均划分为ε段。因此，向量a可表示为:
Figure BDA0002041337110000142
其中：D为铣刀直径；该三维稳定性模型的叶瓣图如图3所示。
所述步骤三中，根据数控机床铣削过程建立以高效、低成本以及高利润为对象的铣削工艺参数优化模型的目标函数包含以下目标函数：
(3-1)建立铣削耗时的目标函数：
Figure BDA0002041337110000143
其中，D为铣刀直径；Ωri为第i次粗加工时机床的主轴转速，Ωs为精加工时机床的主轴转速，n为粗加工的次数，z为铣刀齿数，fri为第i次粗加工时候每齿进给量；fs为精加工时每齿进给量；h1以及h2为与走刀长度以及切入、切出时间相关的常数；Ttc为换刀时间；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命，Lri为第i次粗加工时的走刀长度，LS为精加工时的走刀长度；
Figure BDA0002041337110000151
Figure BDA0002041337110000152
Figure BDA0002041337110000153
Figure BDA0002041337110000154
其中dri为第i次粗加工时机床的轴向切深；ari为第i次粗加工时机床径向切深，Ωs为精加工时机床的主轴转速；as为精加工时机床的径向切；Cv,Kv,xv,yv,sv,qv,pv和l是与刀具和工件材料相关的常数；GInt()和SInt()分别表示取最大整数和最小整数，W和L分别为工件的宽度以及长度，e’是为避免撞刀而设置的任意距离；
(3-2)建立铣削成本的目标函数：
Figure BDA0002041337110000155
其中：k0为单位时间的人工及管理成本，kt为单位时间的加权刀具成本，Tp为加工准备时间。
(3-3)建立铣削利润率的目标函数：
Figure BDA0002041337110000156
其中：Sp为产品的市场价格，Cmat为原材料成本。
所述步骤四中，根据数控机床性能、刀具性能、工件性能、加工要求以及步骤二中的三维稳定性要求建立铣削工艺参数优化模型的约束条件包含以下约束条件：
(3-4)主轴转速约束条件
Ωmin＜Ωri＜Ωmax；Ωmin＜Ωs＜Ωmax；
其中，Ωmax以及Ωmin为机床的最大、最小主轴转速；Ωri为第i次粗加工时机床的主轴转速；Ωs为精加工时机床的主轴转速。
(3-5)进给速度约束条件
fmin＜fri＜fmax；fmin＜fs＜fmax；
其中，fmax以及fmin为机床的最大、最小每齿进给量即进给速度；fri为第i次粗加工时机床的每齿进给量；fs为精加工时机床的每齿进给量；
(3-6)轴向切深约束条件
dmin＜dri＜dmax；dmin＜ds＜dmax；
其中，dmax以及dmin为机床的最大、最小轴向切深；dri为第i次粗加工时机床的轴向切深；ds为精加工时机床的轴向切深；dt为总铣削深度；n为粗加工的总次数。
除此之外，对于多次走刀端铣而言，必须满足：
Figure BDA0002041337110000161
其中：dt为总铣削深度。
(3-7)径向切深约束条件
amin＜ari＜amax；amin＜as＜amax；
其中，amax以及amin为机床的最大、最小径向切深；ari为第i次粗加工时机床的径向切深；as为精加工时机床的径向切深。
(3-8)切削力约束条件
Figure BDA0002041337110000162
其中，Cu,Ku,pu,qu,su,xu以及yu为常数，Fmin以及Fmax为机床的最大、最小主切削力；Fri为第i次粗加工时机床的主切削力；Fs为精加工时机床的主切削力。
(3-9)切削力矩约束条件
Figure BDA0002041337110000163
其中：TMmax为机床的最大主轴转矩；TMri为第i次粗加工时机床的主轴转矩；TMs为精加工时机床的主轴转矩。
(3-10)主轴功率约束条件
Figure BDA0002041337110000164
其中：pmax为机床的最大主轴功率；Pri为第i次粗加工时机床的主轴功率；Ps为精加工时机床的主轴功率。
(3-11)表面粗糙度约束条件
Figure BDA0002041337110000171
其中，re为刀尖圆角半径，Rrmax以及Rsmax分别为粗加工和精加工最大表面粗糙度；Rri为第i次粗加工时表面粗糙度；Rs为精加工时表面粗糙度。
(3-12)刀具寿命约束条件
述刀具寿命约束条件为：
tri≥TR，ts≥TR；
其中，TR为使用者所要求的刀具最短加工寿命；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命。
(3-13)稳定性约束条件
如步骤二所建立的铣削稳定性模型，主轴转速Ω，轴向切深d以及径向切深a必须满足三维颤振稳定性模型所要求的临界条件。
所述步骤五的具体流程如图1所示，
步骤5.1：对每只萤火虫进行编码，并根据所设定编码方案和编码中的上边界L和下边界U初始化每只萤火虫的坐标，具体编码方式如下所示：
pA＝(n，dr，dh，ds，ar，ah，as，fr，fh，fs，Ωr，Ωh，Ωs)
L＝(nmin，drmin，dhmin，dsmin，armin，ahmin，asmin，frmin，fhmin，fsmin，Ωrmin，Ωhmin，Ωsmin)
U＝(nmax，drmax，dhmax，dsmax，armax，ahmax，asmax，frmax，fhmax，fsmax，Ωrmax，Ωhmax，Ωsmax)
其中，pA为第A只萤火虫的编码，n为粗加工的次数，dh为半精加工时的轴向切深，ah为半精加工时径向切深，fh为半径加工时的进给速度；Ωh为半精加工时主轴转速；
步骤:5.2：计算初始种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度G(pA)；所述约束违背程度G(pA)的具体计算方法如下所示：
Figure BDA0002041337110000172
GC(pA)＝max{0，gC(pA)}
其中,G(pA)为第A个萤火虫的总约束违背程度，GC(pA)为第A个萤火虫的第C个约束条件的违背程度；CMAX为约束条件的总数量；若在计算时存在等式约束H(p)，则将其转变为不等式约束：|H(p)|-δ≤0,其中δ＝1.0E-004；
步骤5.3：根据选择Pareto前沿解的标准在当前萤火虫种群中获取Pareto前沿解，将所述Pareto前沿解放入外部存档中，如果外部存档超过预设的范围，则利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤5.4：在当前Pareto前沿解中选择最优萤火虫；具体选择方法为：将Pareto前沿解中的所有萤火虫按照拥挤度的升序排序，并在排序的顶部选择一定比例的萤火虫作为最优萤火虫候选，将候选的萤火虫按照每个目标函数值的大小进行排序，选择任意两个毗邻的萤火虫，以这两个萤火虫坐标连线的中点作为最优萤火虫的坐标；
步骤5.5：比较当前种群中任意两个萤火虫A与H的支配关系，若A支配H，则转步骤5.6；否则转步骤5.7；
步骤5.6：利用基于云模型的坐标摄动机制处理萤火虫A的坐标，利用基于拥挤度的自适应吸引力机制更新萤火虫H的坐标；并转向步骤5.8；
步骤5.7：利用基于云模型的坐标摄动机制处理最优萤火虫的坐标，利用基于拥挤度的最优萤火虫引导机制更新萤火虫H坐标；
步骤5.8：移动过之后，按照一定的比例选取萤火虫，再次利用基于云模型的坐标摄动机制对选取的萤火虫的坐标进行处理；计算当前种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；并转步骤5.9和步骤5.10；
步骤5.9：根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，将其暂时存入外部存档，并转向步骤5.12；
步骤5.10：利用基于云模型的坐标摄动机制对步骤3.8获得的种群中的每只萤火虫的坐标进行处理，并计算处理后的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；
步骤5.11：再次根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，将其暂时存入外部存档；
步骤5.12：对外部存档中的所有的Pareto前沿解按照选择Pareto前沿解的标准，重新获取Pareto前沿解，并剔除非前沿解；
步骤5.13：如果外部存档超过预设的范围，利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤5.14：IterNum＝IterNum+1；IterNum＞IterNummax？如果是转步骤5.15，否则转步骤5.4；
步骤5.15：输出外部存档中的所有萤火虫，即多组具有代表性的Pareto前沿解。
在每次迭代计算之前基于聚类原理多目标萤火虫优化算法会根据总切削深度dt，针对多次走刀高速铣削中的粗加工、半精加工以及精加工的轴向切削的深度进行划分，具体划分方案为：
Figure BDA0002041337110000191
则前n次走刀为粗加工,第n+1次走刀为精加工,如果n不是整数，则粗加工次数为；
Figure BDA0002041337110000192
dh＝dt-ds-n*dr
则前n次走刀为粗加工，第n+1次走刀为半精加工；第n+2次走刀为精加工。本实施例的走刀路线图如图2所示。
所述选择Pareto前沿解的标准为：将每个萤火虫与其他萤火虫依次组成一对，如果两个萤火虫均为可行解即G(pA)＝0，在该对中选择占主导地位的萤火虫为该对的前沿解，如果两者均不占主导地位，则随机选择一个萤火虫作为该对的前沿解，如果一对中只有一个可行解则选择可行解为该对的前沿解，如果两者均不是可行解，则选择G(pA)最小的萤火虫作为该对的前沿解，所述占主导地位的萤火虫为Ttotal值小、Ctotal值小、ptotal值大的萤火虫。
为了保证优化所得结果具有良好的多样性及分散性。本实施例采用了一种经典的拥挤度计算方法来估计特定解周围解的密度。该方法的详细描述如下：
拥挤度计算需要根据对所有萤火虫按照每个目标函数值按升序进行排序。任何萤火虫的拥挤度计算值都可以通过该萤火虫在每个目标函数方向上毗邻个体的平均距离来计算。以Pareto前沿解中的第A个萤火虫为例，如附图4所示。拥挤度计算值CrowdA为矩形的平均边长，即：
Figure BDA0002041337110000193
其中Emax为目标函数总数。fE，A+1以及fE，A-1分别表示第A+1及第A-1个萤火虫的第E个目标函数值。CrowdA越小，则该萤火虫周围就越拥挤，反之则越疏散。具有任意最大或最小目标函数的萤火虫的拥挤度设定为无穷大。
所述基于拥挤度的自适应吸引力机制为：
萤火虫算法的优化机制灵感来自萤火虫的闪光行为，它遵循以下三个原则：(1)萤火虫是雌雄同体的，因此萤火虫会不区分性别的被其他萤火虫所吸引；(2)萤火虫的吸引力与其亮度成正比且随着距离的增加而减小，因此，对于任何两只闪烁的萤火虫，较暗的那只就会被吸引且飞向较亮的那只，如果没有萤火虫比该萤火虫更亮，则它将会随机移动；(3)萤火虫的亮度由其目标函数的大小决定。在上述三个原则的基础上，将拥挤度引入吸引力因子中，进一步提高该算法得到的Pareto前沿解的多样性。以萤火虫A和H为例，假定H相对较暗而被A所吸引，具体操作如下:
Figure BDA0002041337110000201
其中，
Figure BDA0002041337110000202
为本次迭代计算中萤火虫H被基于拥挤度的自适应吸引力机制处理后的坐标；
Figure BDA0002041337110000203
本次迭代计算中萤火虫H被处理前的坐标；
Figure BDA0002041337110000204
为萤火虫A的坐标；
Figure BDA0002041337110000205
为萤火A对萤火H的吸引度，β0为r＝0时萤火虫的吸引力；γ为光吸收系数；α′为随机化参数；
Figure BDA0002041337110000206
为与
Figure BDA0002041337110000207
大小成正相关的随机数向量；Rand(SizeVar)为[0,1]之间的SizeVar大小的随机数矩阵，SizeVar设计变量的规模，Distance(·)表示求解的距离；CrowdH为萤火虫H的拥挤度；
通过将目标萤火虫的拥挤度引入到吸引力因子中，萤火虫则更倾向于飞向密度较小的区域，即当目标萤火虫周围的密度较大时，则吸引力相对较弱，则萤火虫飞行动作中随机性所占的比例相对较大，相反，如果目标萤火虫周围的密度较小，则吸引力较强，那么萤火虫更倾向于飞向目标萤火虫。
所述基于拥挤度的最优萤火虫引导机制如下所示：
优化进行到一定阶段后，很有可能出现Pareto前沿解互不支配的情况。在这种情况下，在传统萤火虫算法中的萤火虫会随机移动，这种随机优化机制效率极低，甚至会导致一些萤火虫陷入毫无意义的随机徘徊而永不收敛。因此，本专利还提出了一种基于拥挤度的最优萤火虫引导机制。在这个机制下，如果Pareto前沿解中所有萤火虫均不互相支配，则他们将会向规定的最优位置飞行。以萤火虫H为例，具体操作如下。
Figure BDA0002041337110000211
其中，
Figure BDA0002041337110000212
为第IterNum次迭代计算中萤火H被基于拥挤度的最优萤火虫引导机制处理后的坐标，
Figure BDA0002041337110000213
为萤火H被处理前的坐标，
Figure BDA0002041337110000214
为最优萤火虫对萤火虫H的吸引力，
Figure BDA0002041337110000215
为最优萤火虫的坐标；
为了避免所提出的萤火虫算法在其优化过程中可能出现的过早收敛的问题，基于正态云模型可在定性概念与定量值之间不确定性转换的优良特性，提出了基于云模型的坐标摄动机制。该机制的详细描述如下:
Figure BDA0002041337110000216
其中，En以及He分别表示随着迭代次数的增加而减小的熵和超熵，
Figure BDA0002041337110000217
表示被基于云模型的坐标摄动机制处理前萤火中A的坐标，
Figure BDA0002041337110000218
表示被处理后的萤虫A的坐标；
Figure BDA0002041337110000219
表示被基于云模型的坐标摄动机制处理前最优萤火中O的坐标；
Figure BDA00020413371100002110
表示被处理后最优萤火中O的坐标。
基于聚类原理的外部存档更新机制
开辟外部存档保存搜索到的相对较为优秀的Pareto前沿解，将每代中的Pareto前沿解与外部存档的内容进行对比，以收集当代中的Pareto前沿解，并剔除存储在存档中的劣解，且删除所有被新存入存档中的解所主导的所有解。一般来说，外部存储是根据支配程度和拥挤度进行更新的，这可能会导致局部过度拥挤的情况下多样性减小的问题。为了避免这些缺陷，本专利提出了一种基于聚类原理的外部存档更新策略。本专利中外部存档中萤火虫均为Pareto前沿解，所以更新的主要目的是为了保证种群的多样性。所提出的机制不仅能在保证多样性的前提下有效地减少外部存档，而且能保证所选萤火虫具有一定的代表性，该方法的实现过程总结如下:
步骤S1：初始化集群集Q：外部存档中的每个Pareto前沿解组成一个集群。
步骤S2：WHILE集群集中的集群数量没有减少到预先指定的数量
DO
S1.1.计算所有集群对之间的距离：两个集群Cδ∈Q，Cβ∈Q之间的距离为两个集群中萤火虫对之间的平均距离:
Figure BDA0002041337110000221
其中||·||表示萤火虫集群之间的距离。
S1.2.确定距离最小的两个集群Cδ以及Cβ，然后以这两个集群的并集置替换为这两个集群Q＝Q-{Cδ，Cβ}∪{Cδ∪Cβ}。
步骤S3：END WHILE
步骤S4：FOR集群集C中的每个集群
IF集群中萤火虫的数量为1
选择这仅有的萤火虫作为代表。
ELSEIF集群中萤火虫数量为2
选择拥挤度较大的萤火虫为代表。
ELSE IF集群中萤火虫数量超过2
选择与集群中其他萤火虫平均距离最小的萤火虫作为代表。
步骤S5：END IF
步骤S6：END FOR
步骤S7：输出消减后的外部存档。
所述伪权重估计法为：
在找到多个具有代表性的解后，为了方便决策者选择一个具有代表性的解，采用了伪权系数向量法对每个典型解计算了一个确定的不同目标重要性的伪权系数向量。对于目标函数的最小化，对于得到的集合中的任何解，每个目标函数i的权值wi计算如下:
Figure BDA0002041337110000231
其中，wy为某个前沿解的第y个目标函数的权重值，fy,max为所有前沿解的第y个函数的最大值；fy为该某个前沿解第y个目标函数的值；fy,min为所有前沿解的第y个函数的最小值；Emax为所有目标函数的总个数，fE,max为所有前沿解的第E个函数的最大值，fE为该某个前沿解的第E个目标函数的值，fE,min为所有前沿解的第E个函数的最小值。若一个有代表性的解接近个体目标函数的最小值，则这个函数的权重值很大，设置为1。对于目标函数为最大值的情况，(fy，max-fy)需要替代为(fy-fy，min)。
如图5所示，为本发明所提出的算法进行多目标多约束的铣削参数优化的求解结果以及决策结果，图(a)为优化算法求解所得关于加工时长、加工成本以及利润率三个目标的500组Pareto前沿解；图(b)为通过聚类算法决策更新消减之后所划分出的20组种群，由图可看出更新消减后所得的20组解分布较为均匀；由目标函数值的伪权重算法可看出20组解的代表性较强，具体20组解的参数如表1所示；
表1
Figure BDA0002041337110000232
Figure BDA0002041337110000241
Figure BDA0002041337110000251
如图6所示，为本实施例所提出的算法在六个标准测试函数CONSTR、TNK、KITA、SRN、BNH以及OSY上的优化结果，可见所提出的算法可有效的收敛到Pareto最优解，且均匀性相对较好。表2为本实施例所提出算法在六个标准测试函数：CONSTR、TNK、KITA、SRN、BNH以及OSY上的优化程度、均匀性和延展性。
表2
Figure BDA0002041337110000252
Figure BDA0002041337110000261

Claims (7)

1.一种引入三维稳定性约束的铣削加工多目标优化决策方法，其特征在于，包括如下步骤：
步骤1：基于二阶全离散法，在主轴转速Ω、轴向切深d的铣削稳定性预测模型中引入径向切深a，构建关于主轴转速Ω、轴向切深d以及径向切深a的铣削三维稳定性预测模型；
步骤2：建立多目标多约束的多次走刀铣削参数优化模型，该模型包括多目标函数；和根据数控机床性能、刀具性能、工件性能、加工要求以及步骤1中的三维稳定性模型建立的多约束条件；
步骤3：利用基于聚类原理的多目标萤火虫优化算法求解多目标多约束的多次走刀铣削参数优化模型，得到多组具有代表性的切削参数，即多组具有代表性的Pareto前沿解，所述基于聚类原理的多目标萤火虫优化算法包括：基于云模型的坐标摄动机制、基于拥挤度的自适应吸引力机制、基于拥挤度的最优萤火虫引导机制、基于聚类原理的外部存档更新消减机制；
步骤4：利用伪权重估计法，计算每个具有代表性的前沿解其不同目标函数的权重值；
所述步骤2中的多目标函数包括：铣削耗时的目标函数、铣削成本的目标函数、铣削利润的目标函数；
所述铣削耗时的目标函数如下所示：
Figure RE-FDA0003712620150000011
其中，D为铣刀直径；Ωri为第i次粗加工时机床的主轴转速，Ωs为精加工时机床的主轴转速，n为粗加工的次数，z为铣刀齿数，fri为第i次粗加工时候每齿进给量；fs为精加工时每齿进给量；h1以及h2为与走刀长度以及切入、切出时间相关的常数；Ttc为换刀时间；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命，Lri为第i次粗加工时的走刀长度，LS为精加工时的走刀长度；
Figure RE-FDA0003712620150000012
Figure RE-FDA0003712620150000021
Figure RE-FDA0003712620150000022
Figure RE-FDA0003712620150000023
其中dri为第i次粗加工时机床的轴向切深；ari为第i次粗加工时机床径向切深，Ωs为精加工时机床的主轴转速；as为精加工时机床的径向切；Cv,Kv,xv,yv,sv,qv,pv和l是与刀具和工件材料相关的常数；GInt( )和SInt( )分别表示取最大整数和最小整数，W和L分别为工件的宽度以及长度，e’是为避免撞刀而设置的任意距离；
所述铣削成本的目标函数如下所示：
Figure RE-FDA0003712620150000024
其中：k0为单位时间的人工及管理成本，kt为单位时间的加权刀具成本，Tp为加工准备时间；
所述铣削利润的目标函数如下所示：
Figure RE-FDA0003712620150000025
其中：Sp为产品的市场价格，Cmat为原材料成本；
所述多约束条件包括：主轴转速约束条件、进给速度约束条件、轴向切深约束条件、径向切深约束条件、切削力约束条件、切削力矩约束条件、主轴功率约束条件、表面粗糙度约束条件、刀具寿命约束条件和三维稳定性模型构成的约束条件；
所述主轴转速约束条件为：
Ωmin＜Ωri＜Ωmax；Ωmin＜Ωs＜Ωmax；
其中，Ωmax以及Ωmin为机床的最大、最小主轴转速；Ωri为第i次粗加工时机床的主轴转速；Ωs为精加工时机床的主轴转速；
所述进给速度约束条件为：
fmin＜fri＜fmax；fmin＜fs＜fmax；
其中，fmax以及fmin为机床的最大、最小每齿进给量即进给速度；fri为第i次粗加工时机床的每齿进给量；fs为精加工时机床的每齿进给量；
所述轴向切深约束条件为：
dmin＜dri＜dmax；dmin＜ds＜dmax；
Figure RE-FDA0003712620150000031
其中，dmax以及dmin为机床的最大、最小轴向切深；dri为第i次粗加工时机床的轴向切深；ds为精加工时机床的轴向切深；dt为总铣削深度；n为粗加工的总次数；
所述径向切深约束条件为：
amin＜ari＜amax；amin＜as＜amax；
其中，amax以及amin为机床的最大、最小径向切深；ari为第i次粗加工时机床的径向切深；as为精加工时机床的径向切深；
所示切削力约束条件为：
Figure RE-FDA0003712620150000032
其中，Cu，Ku，pu，qu，su，xu以及yu为常数，Fmin以及Fmax为机床的最大、最小主切削力；Fri为第i次粗加工时机床的主切削力；Fs为精加工时机床的主切削力；
所述切削力矩约束条件为：
Figure RE-FDA0003712620150000033
其中：TMmax为机床的最大主轴转矩；TMri为第i次粗加工时机床的主轴转矩；TMs为精加工时机床的主轴转矩；
所述主轴功率约束条件为：
Figure RE-FDA0003712620150000034
其中：pmax为机床的最大主轴功率；Pri为第i次粗加工时机床的主轴功率；Ps为精加工时机床的主轴功率；
所述表面粗糙度约束条件为：
Figure RE-FDA0003712620150000041
其中，re为刀尖圆角半径，Rrmax以及Rsmax分别为粗加工和精加工最大表面粗糙度；Rri为第i次粗加工时表面粗糙度；Rs为精加工时表面粗糙度；
所述刀具寿命约束条件为：
tri≥TR，ts≥TR；
其中，TR为使用者所要求的刀具最短加工寿命；tri为第i次粗加工时刀具的寿命；ts为精加工时刀具的寿命；
所述三维稳定性模型构成的约束条件为求解步骤二中所建立的铣削稳定性模型时，主轴转速Ω，轴向切深d以及径向切深a满足三维颤振稳定性模型所要求的临界条件；
所述得到若干个具有代表性的解的具体步骤如下所示：
步骤3.1：对每只萤火虫进行编码，并根据所设定编码方案和编码中的上边界L和下边界U初始化每只萤火虫的坐标，具体编码方式如下所示：
pA＝(n，dr，dh，ds，ar，ah，as，fr，fh，fs，Ωr，Ωh，Ωs)
L＝(nmin，drmin，dhmin，dsmin，armin，ahmin，asmin，frmin，fhmin，fsmin，Ωrmin，Ωhmin，Ωsmin)
U＝(nmax，drmax，dhmax，dsmax，armax，ahmax，asmax，frmax，fhmax，fsmax，Ωrmax，Ωhmax，Ωsmax)
其中，pA为第A只萤火虫的编码，n为粗加工的次数，dh为半精加工时的轴向切深，ah为半精加工时径向切深，fh为半精加工时的进给速度；Ωh为半精加工时主轴转速；
步骤3.2：计算初始种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度G(pA)；所述约束违背程度G(pA)的具体计算方法如下所示：
Figure RE-FDA0003712620150000042
GC(pA)＝max{0，gC(pA)}
其中，G(pA)为第A个萤火虫的总约束违背程度，GC(pA)为第A个萤火虫的第C个约束条件的违背程度；CMAX为约束条件的总数量；若在计算时存在等式约束H(p)，则将其转变为不等式约束：|H(p)|-δ≤0，其中δ＝1.0E-004；
步骤3.3：根据选择Pareto前沿解的标准在当前萤火虫种群中获取Pareto前沿解，将所述Pareto前沿解放入外部存档中，如果外部存档超过预设的范围，则利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤3.4：在当前Pareto前沿解中选择最优萤火虫；具体选择方法为：将Pareto前沿解中的所有萤火虫按照拥挤度的升序排序，并在排序的顶部选择一定比例的萤火虫作为最优萤火虫候选，将候选的萤火虫按照每个目标函数值的大小进行排序，选择任意两个毗邻的萤火虫，以这两个萤火虫坐标连线的中点作为最优萤火虫的坐标；
步骤3.5：比较当前种群中任意两个萤火虫A与H的支配关系，若A支配H，则转步骤3.6；否则转步骤3.7；
步骤3.6：利用基于云模型的坐标摄动机制处理萤火虫A的坐标，利用基于拥挤度的自适应吸引力机制更新萤火虫H的坐标；并转向步骤3.8；
步骤3.7：利用基于云模型的坐标摄动机制处理最优萤火虫的坐标，利用基于拥挤度的最优萤火虫引导机制更新萤火虫H坐标；
步骤3.8：每个萤火虫的坐标均被移动过之后，按照一定的比例选取萤火虫，再次利用基于云模型的坐标摄动机制对选取的萤火虫的坐标进行处理；计算当前种群中的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；并转步骤3.9和步骤3.10；
步骤3.9：根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，将其暂时存入外部存档，并转向步骤3.12；
步骤3.10：利用基于云模型的坐标摄动机制对步骤3.8获得的种群中的每只萤火虫的坐标进行处理，并计算处理后的每只萤火虫的多目标函数值、拥挤程度以及约束违背程度；
步骤3.11：再次根据选择Pareto前沿解的标准获取当前种群中的Pareto前沿解，并将其暂时存入外部存档；
步骤3.12：对外部存档中的所有的Pareto前沿解按照选择Pareto前沿解的标准，重新获取Pareto前沿解，并剔除非前沿解；
步骤3.13：如果外部存档超过预设的范围，利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减；
步骤3.14：IterNum＝IterNum+1；判断IterNum是否大于IterNummax，如果是转步骤3.15，否则转步骤3.4；
步骤3.15：输出外部存档中的所有萤火虫，即多组具有代表性的Pareto前沿解。
2.根据权利要求1所述的方法，其特征在于，所述径向切深a的表达式为：
Figure RE-FDA0003712620150000061
其中s为无穷小量，ε为径向切深被平均划分为ε段，D为铣刀直径。
3.根据权利要求1所述的方法，其特征在于，所述计算具有代表性的前沿解其不同目标函数的权重值的具体方法如下所示：
Figure RE-FDA0003712620150000062
其中，wy为某个前沿解的第y个目标函数的权重值，fy,max为所有前沿解的第y个函数的最大值；fy为该某个前沿解第y个目标函数的值；fy,min为所有前沿解的第y个函数的最小值；Emax为所有目标函数的总个数，fE,max为所有前沿解的第E个函数的最大值，fE为该某个前沿解的第E个目标函数的值，fE,min为所有前沿解的第E个函数的最小值。
4.根据权利要求1所述的方法，其特征在于，在每次迭代计算之前基于聚类原理的多目标萤火虫优化算法会根据总切削深度dt，针对多次走刀高速铣削中的粗加工、半精加工以及精加工的轴向切削的深度进行划分，具体划分方案为：
Figure RE-FDA0003712620150000063
则前n次走刀为粗加工,第n+1次走刀为精加工,如果n不是整数，则粗加工次数为；
Figure RE-FDA0003712620150000064
dh＝dt-ds-n*dr
则前n次走刀为粗加工，第n+1次走刀为半精加工；第n+2次走刀为精加工。
5.根据权利要求1所述的方法，其特征在于，所述选择Pareto前沿解的标准为：将每个萤火虫与其他萤火虫依次组成一对，如果两个萤火虫均为可行解即G(pA)＝0，在该对中选择占主导地位的萤火虫为该对的前沿解，如果两者均不占主导地位，则随机选择一个萤火虫作为该对的前沿解，如果一对中只有一个可行解则选择可行解为该对的前沿解，如果两者均不是可行解，则选择G(pA)最小的萤火虫作为该对的前沿解，所述占主导地位的萤火虫为Ttotal值小、Ctotal值小、ptotal值大的萤火虫。
6.根据权利要求1所述的方法，其特征在于，所述基于拥挤度的自适应吸引力机制为：
Figure RE-FDA0003712620150000071
其中，
Figure RE-FDA0003712620150000072
为本次迭代计算中萤火虫H被基于拥挤度的自适应吸引力机制处理后的坐标；
Figure RE-FDA0003712620150000073
本次迭代计算中萤火虫H被处理前的坐标；
Figure RE-FDA0003712620150000074
为萤火虫A的坐标；
Figure RE-FDA0003712620150000075
为萤火A对萤火H的吸引度，β0为r＝0时萤火虫的吸引力；γ为光吸收系数；α'为随机化参数；
Figure RE-FDA0003712620150000076
为与
Figure RE-FDA0003712620150000077
大小成正相关的随机数向量；Rand(SizeVar)为[0,1]之间的SizeVar大小的随机数矩阵，SizeVar为设计变量的规模，Distance(·)表示求解的距离；CrowdH为萤火虫H的拥挤度；
所述基于拥挤度的最优萤火虫引导机制如下所示：
Figure RE-FDA0003712620150000078
其中，
Figure RE-FDA0003712620150000079
为第IterNum次迭代计算中萤火H被基于拥挤度的最优萤火虫引导机制处理后的坐标，
Figure RE-FDA00037126201500000710
为萤火H被处理前的坐标，
Figure RE-FDA00037126201500000711
为最优萤火虫对萤火虫H的吸引力，
Figure RE-FDA00037126201500000712
为最优萤火虫的坐标；
所述基于云模型的坐标摄动机制的公式如下所示：
Figure RE-FDA0003712620150000081
其中，En以及He分别表示随着迭代次数的增加而减小的熵和超熵，
Figure RE-FDA0003712620150000082
表示被基于云模型的坐标摄动机制处理前萤火中A的坐标，
Figure RE-FDA0003712620150000083
表示被处理后的萤虫A的坐标；
Figure RE-FDA0003712620150000084
表示被基于云模型的坐标摄动机制处理前最优萤火O的坐标；
Figure RE-FDA0003712620150000085
表示被处理后最优萤火O的坐标。
7.根据权利要求1所述的方法，其特征在于,所述利用基于聚类原理的外部存档更新消减机制对外部存档进行更新消减的具体操作如下所示：
步骤A：初始化集群集Q，该集群集中包括本次迭代计算得到的所有的Pareto前沿解；每个Pareto前沿解构成一个独立的集群；
步骤B：计算集群集中，每个集群与其他集群之间的距离d的计算公式如下所示，
Figure RE-FDA0003712620150000086
其中：||·||代表为欧式距离；I1为集群Cδ中的萤火虫，I2为集群Cβ中的萤火虫，Cδ∈Q,Cβ∈Q；
步骤C：选择两个距离最小的集群，并用这两个集群的并集代替这两个集群，则集群集Q的变化如下所示：
Q＝Q-{Cδ，Cβ}∪{Cδ∪Cβ}
步骤D：并计算该集群集Q中集群的数量是否减少到预先指定的数量，如果是则转步骤E；否则转步骤B；
步骤E：最终得到指定数量的集群，判断每个集群中的萤火虫数量是否为1，如果集群中萤火虫的数量为2，则选择拥挤度大的萤火为这个集群的代表，如果数量超过2 选择与该集群中其他萤火虫平均距离最小的萤火虫作为该集群的代表。
CN201910342927.6A 2019-04-26 2019-04-26 一种引入三维稳定性约束的铣削加工多目标优化决策方法 Active CN110162841B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910342927.6A CN110162841B (zh) 2019-04-26 2019-04-26 一种引入三维稳定性约束的铣削加工多目标优化决策方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910342927.6A CN110162841B (zh) 2019-04-26 2019-04-26 一种引入三维稳定性约束的铣削加工多目标优化决策方法

Publications (2)

Publication Number Publication Date
CN110162841A CN110162841A (zh) 2019-08-23
CN110162841B true CN110162841B (zh) 2022-09-13

Family

ID=67640039

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910342927.6A Active CN110162841B (zh) 2019-04-26 2019-04-26 一种引入三维稳定性约束的铣削加工多目标优化决策方法

Country Status (1)

Country Link
CN (1) CN110162841B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112296357B (zh) * 2020-10-29 2021-11-30 福州大学 面向激光选区熔融工艺的增材制造相同并行机调度方法
CN114327859B (zh) * 2021-11-18 2024-06-07 西安电子科技大学 云计算环境大规模问题代理优化的源模型聚类选择方法
CN114939772B (zh) * 2022-04-20 2023-12-19 湖南星途航空航天器制造有限公司 一种加工薄壁弧形件的力适应性方法
CN116662859B (zh) * 2023-05-31 2024-04-19 西安工程大学 非遗文化数据特征选择方法
CN117389161B (zh) * 2023-12-12 2024-02-27 山东科技大学 考虑多加工阶段的平面立铣参数优化及加工性能预测方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102609591A (zh) * 2012-02-16 2012-07-25 华中科技大学 一种重型机床切削参数的优化方法
CN106257849A (zh) * 2016-09-13 2016-12-28 哈尔滨工程大学 基于多目标量子萤火虫搜索机制的频谱感知方法
CN108416474A (zh) * 2018-03-02 2018-08-17 南京工程学院 基于虚拟适应度和Pareto解集分层的改进优化算法
CN108828436A (zh) * 2018-06-27 2018-11-16 桂林电子科技大学 基于混沌云自适应萤火虫算法的模拟电路故障诊断方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012107594A1 (en) * 2011-02-11 2012-08-16 Ecole Polytechnique Federale De Lausanne (Epfl) High speed pocket milling optimisation

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102609591A (zh) * 2012-02-16 2012-07-25 华中科技大学 一种重型机床切削参数的优化方法
CN106257849A (zh) * 2016-09-13 2016-12-28 哈尔滨工程大学 基于多目标量子萤火虫搜索机制的频谱感知方法
CN108416474A (zh) * 2018-03-02 2018-08-17 南京工程学院 基于虚拟适应度和Pareto解集分层的改进优化算法
CN108828436A (zh) * 2018-06-27 2018-11-16 桂林电子科技大学 基于混沌云自适应萤火虫算法的模拟电路故障诊断方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Goran R Miodragovic＇等.Optimization of multi-pass turning and multi-pass face milling using subpopulation firefly algorithm.《institution of mechanical engineers》.2018, *
一种面向绿色高效的数控铣削参数优化方法;王培建等;《航空制造技术》;20160401(第07期);第50-54页 *
考虑线缆干涉约束的铣削机器人进刀轨迹规划;罗龙健;《中国优秀博硕士学位论文全文数据库(硕士) 信息科技辑》;20180215(第02期);全文 *
高速铣削刀轨优化技术的研究;孙全平 等;《机械科学与技术》;20040831;第23卷(第8期);第922-926页 *

Also Published As

Publication number Publication date
CN110162841A (zh) 2019-08-23

Similar Documents

Publication Publication Date Title
CN110162841B (zh) 一种引入三维稳定性约束的铣削加工多目标优化决策方法
CN107480352B (zh) 一种铣削加工工艺参数的可靠性优化方法
CN111563301A (zh) 一种薄壁件铣削加工参数优化方法
CN106485339A (zh) 一种电力系统的负荷特性确定方法及系统
CN112380760B (zh) 基于多算法融合的多目标工艺参数智能优化方法
Yang et al. Multi-objective optimization of multi-pass face milling using particle swarm intelligence
CN113050543B (zh) 一种切削阶段数控车床加工参数节能优化方法
CN104808629A (zh) 一种柔性机器人制造单元的调度方法
CN106774162A (zh) 一种数控加工参数多目标优化方法
Bhuiyan et al. Optimization of cutting parameters in turning process
CN114415595B (zh) 一种车削优化方法、系统、计算机设备和存储介质
CN113960971B (zh) 一种基于行为决策网络粒子群优化的柔性车间调度方法
CN108829036B (zh) 一种金属铸坯切削成形加工过程的优化调度方法
Pytlak Multicriteria optimization of hard turning operation of the hardened 18HGT steel
Bhuiyan et al. Optimization of cutting parameters in turning process
Shin et al. A comparative study of PSO, GSA and SCA in parameters optimization of surface grinding process
Osyczka et al. Optimum design of robot grippers using genetic algorithms
CN108876654B (zh) 一种多类电缆加工的优化调度方法
Park et al. Multi-objective optimization of turning process for hardened material based on hybrid approach
Wang et al. Multi-objective optimization of high-speed milling with parallel genetic simulated annealing
CN114734210B (zh) 一种基于滚动填炉策略的钢管合同组批方法
CN115808909A (zh) 一种面向模具热处理准时与节能生产的动态批调度方法
Xi et al. Cutting parameter optimization based on particle swarm optimization
CN114237166A (zh) 基于改进spea2算法求解多转速节能调度问题的方法
Kumar et al. Integrated total cost and tolerance optimization with genetic algorithm

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant