基于牛顿迭代算法的伺服压力机控制系统研究.pdf

第41卷第12期Vol． 41 No． 12FORGING 利用MATLAB软件拟合目标工艺曲线,基于牛顿迭代算法实现了利用滑块位移对曲柄转角的多点一次性精确求解,并进一步利用计算结果拟合出曲柄转角-时间的函数;推导出控制周期内目标工艺曲线滑块位移对应的脉冲频率与个数,利用STM32实现输出目标脉冲。研究结果表明,目标工艺曲线与滑块实际位移曲线基本吻合,证明了该控制系统的设计方案可以实现对滑块位移的精确控制。关键词:伺服压力机;拉深工艺;牛顿迭代算法;运动学方程;滑块位移;曲柄转角DOI: 10． 13330/j． issn． 1000-3940． 2016． 12． 014中图分类号: TP273 文献标识码: A 文章编号: 1000-3940 (2016) 12-0082-06Study on servo press control system based on Newton iteration methodHong Rui, Ye Chunsheng(State Key Laboratory of Materials Processing and Die drawing process; Newton iteration algorithm; kinematical equations; slider block displacement; crank angle收稿日期: 2016 -05 -06;修订日期: 2016 -08 -10基金项目:国家自然科学基金资助项目(51375190)作者简介:洪 瑞(1990 - ),男,硕士E-mail: hongrui_hust@ foxmail． com通讯作者:叶春生(1962 - ),男,博士,副教授E-mail: csyehust@126． com伺服压力机控制系统可以实现对滑块位移的精确控制,目前许多学者已对伺服压力机控制系统进行了深入的研究。文献[1] - [2]介绍了伺服压力机的研究动向,在较早的伺服控制系统设计中,通常采用PC +运动控制板卡、 PC + PLC等模式,之后逐渐发展为采用单片机、 DSP等。伺服压力机包含伺服电机驱动与机械传动,其精确的数学模型计算一直集中于运动机构的逆解运算或永磁同步电机PMSM的控制策略研究。周洋等[3]对四连杆伺服压力机进行运动学逆解分析,并规划滑块的运动轨迹。雷培等[4]建立了四连杆机构的数学模型,利用MATLAB的强大计算和图解功能对其进行运动分析。俞伟洋等[5]求解出在给定滑块工艺曲线下交流永磁同步电机控制模型,得出了电机电枢电流、转速等在周期内随时间变化曲线。以往的研究内容均集中在利用MATLAB进行仿真分析,在利用技术手段实现在复杂工艺曲线条件下,对滑块非线性负载的位万方数据移精确控制方面尚有不足。本文以四连杆压力机为研究对象,以如何实现拉深工艺曲线为例,提出了基于牛顿迭代法伺服控制系统设计方案。1 系统硬件组成伺服控制系统硬件如图1所示,由STM32控制板、接口转换电路、伺服驱动器、伺服电机、开关电源、强电控制电路等组成。 220 V市电经过空气开关QF进入到开关电源和强电控制电路,开关电源完成从交流220 V到直流24 V和5 V电压转换。强电控制电路由交流接触器KM和直流接触器KA,以及一个常闭开关和一个常开开关组成,实现低压24 V对交流220 V的控制。 STM32通过输出方向信号和脉冲信号控制伺服电机,并读取电机速度值和滑块位移,打印输出滑块位移。图1 系统硬件结构Fig． 1 Structure of system hardware2 四连杆运动学分析四连杆结构如图2所示, A点为转轴;曲柄AB是驱动机构,长度为R;连杆BC, CD和CE长度均为L; C点在B点转到水平左侧时达到极限位置,此极限位置处于线段DE上; θ1为曲柄转动角度,逆时针旋转为正; θ2为连杆DC与水平方向夹角,当BC在水平线下方时为正; θ为连杆与垂直方向夹角。在图2中,以A点为原点建立x - y坐标系, E图2 四连杆机构简图Fig． 2 Schematic diagram of four-bar linkage structure点坐标xE和yE为:xE = L - RyE = - Lcosθ + (2L - Lcosθ) = 2L - 2Lcosθ{ (1)求出C点的坐标xC和yC为:xC = L - R + Lsinθ = Rcosθ1 + Lcosθ2yC = L - Lcosθ = Rsinθ1 + Lsinθ2{(2)整理式(1)和式(2)得:2L - 2Lcosθ - yE = 0L - R + Lsinθ - Rcosθ1 - Lcosθ2 = 0L - Lcosθ - Rsinθ1 + Lsinθ2 = 0ìîí (3)将方程组式(3)进行一次求导并整理成矩阵形式:Lcosθ Lsinθ2 0Lsinθ Lcosθ2 02Lsinθ 0 - 1éëùûθ·θ· 2y· Eéëùû=- R θ· 1sinθ1R θ· 1cosθ10éëùû(4)在MATLAB软件中,利用Simulink进行仿真,并根据式(4)编写的矩阵,将曲柄角速度θ· 1以4π· s -1的恒定值作为输入,对该四连杆结构进行仿真。将1个周期(滑块从最高点经过下死点再回到最高点)的仿真结果s, θ, θ2, θ1数据导出,其中s即为E点纵坐标yE,为了方便对比显示,把θ1值缩小10倍,绘制的曲线如图3所示。可以看到, s(t), θ(t)曲线均为抛物线形式,θ2(t)为正弦曲线形式。 s(t)曲线在一个周期内不是单调曲线,考虑将整个周期分为下压阶段(从最高点到下死点)和提升阶段(从下死点回到最高点)两个单调区间,分段进行拟合,拟合出的函数关系38第12期洪 瑞等:基于牛顿迭代算法的伺服压力机控制系统研究 万方数据图3 s(t),θ(t),θ2(t),θ1(t) /10曲线Fig． 3 Curves of s(t),θ(t),θ2(t),θ1(t) /10式如下。(1)下压阶段:θ1(s) = (3 ×10-6) × s6 - 0．0004 × s5 + 0．0203 × s4 -0． 4783 × s3 + 5． 5848 × s2 - 31． 787 × s + 157． 87(5)θ2(s) = - 1 × 10-6 × s6 + 0．0002 × s5 - 0．0075 × s4 +0． 1681 × s3 - 1． 8558 × s2 + 8． 5034 × s + 7． 6137(6)θ(s) = - 8 × 10-5 × s4 + 0． 0071 × s3 -0． 2253 × s2 + 3． 725 × s + 1． 7502 (7)(2)提升阶段:θ1(s) = - 3 × 10-6 × s6 + 0．0004 × s5 - 0．0202 × s4 +0． 4734 × s3 - 5． 5185 × s2 + 32． 341 × s + 202． 44(8)θ2(s) = 1 × 10-6 × s6 - 0．0002 × s5 + 0．0072 × s4 -0． 1627 × s3 + 1． 8243 × s2 - 9． 3132 × s - 7． 7454(9)θ(s) = - 8 × 10-5 × s4 + 0． 007 × s3 - 0． 2215 ×s2 + 3． 6832 × s + 1． 8086 (10)3 基于牛顿迭代法的拉伸工艺下曲柄角度求解本文在采用MATLAB软件曲线拟合工具箱时,为了尽可能减小拟合误差,选用拟合工具箱最高阶拟合,拟合出的拉深工艺曲线下滑块位移傅立叶级数函数为[5 -8]:s(t) = a0 + ∑8i =1(aicosiwt + bisiniwt) (11)式中: a0为第0阶系数; ai, bi分别为第i阶余弦和正弦部分系数; i为阶数; w为角频率,其值设定在6． 284; t为时间。 ai, bi取值如表1所示。根据拟合之后的函数绘制出拉深工艺下滑块位移-时间曲线如图4所示。表1 ai, bi取值Table 1 Parameter values of ai and bii ai bi0 13． 781 13． 13 4． 6812 5． 302 -0． 42643 2． 913 1． 8324 2． 012 0． 70085 1． 068 0． 21326 0． 5584 -0． 11547 0． 2954 -0． 24868 0． 0758 -0． 2893图4 滑块位移-时间曲线Fig． 4 Curve of slider displacement and time采用牛顿迭代法求解非线性方程式(3),基于MATLAB软件利用牛顿迭代法解非线性方程[9 -10]步骤如下。Step1:定义向量x,其中符号变量x1, x2, x3分别代表θ1, θ2, θ,定义非线性方程组矩阵f =[f1 f2 f3],其中: f1 = L - R + L × sinx3 - R × cosx1 - L × cosx2; f2 = L - L × cosx3 - R × sinx1 + L ×sinx2; f3 =2 × L -2 × L × cosx3 - s。Step2:给定初值x0,根的容许误差为e。Step3:计算xk +1 = xk - f′(xk) -1 f(xk),其中,f′(x) = [diff(f,′x1′); diff (f,′x2′); diff (f,′x3′)]Step4:判断x - x0 = e是否成立,如成立转到Step5,否则转到Step3。48锻 压 技 术 第41卷万方数据Step5:迭代结果为x = xk。当t =0． 62 s时,滑块到达下死点,通过分段求解θ1(t)数据点。在下压阶段,即t取0 ~ 0． 62 s时,将t值带入公式(11),求出s值,再将s值带入公式(5), (6), (7),将求得的θ1, θ2, θ值作为初值x0;然后利用基于牛顿迭代法的求解算法解出时间t对应的θ1精确解。同理,在提升阶段,即t取0． 62 ~ 1 s时,将t值带入公式(11),求出s值,再将s值带入公式(8), (9), (10),将求得的θ1, θ2, θ值作为初值x0。通过对滑块位移s(t)曲线以Δt为时间长度进行离散化, s(t)曲线上取200个位移离散点,即Δt =5 ms, s(Δt), s(2Δt), … ,s(nΔt), n =1, 2, … , 200,利用上述牛顿迭代法求得对应200个精确曲柄角度离散点,最终求出曲柄角度-时间曲线,如图5所示。图5 曲柄角度-时间曲线Fig． 5 Curve of crank angle and time利用MATLAB曲线拟合工具箱对得到的θ1(t)数据进行拟合,最终选取偏差最小的七阶多项式拟合曲线,拟合出曲线函数为:θ1(t) = p1 × t7 + p2 × t6 + p3 × t5 + p4 × t4 +p5 × t3 + p6 × t2 + p7 × t + p8 (12)式中: p1 ~ p8为系数,取值如表2所示。表2 多项式系数Table 2 Parameter values of multinomial coefficientp1 p2 p3 p4 p5 p6 p7 p82955 -5745 4152 -4560 5717 -2837 688 3． 127t = nΔt时刻的曲柄角度θ1(nΔt), t = (n + 1)Δt时刻的曲柄角度θ1((n +1)Δt), n = 0,1,2,… ,199,t = nΔt时刻至t = (n + 1)Δt时刻曲柄转速wn +1为:wn+1 = θ1((n + 1)Δt) - θ1(nΔt)Δt (13)同理,得到t = nΔt时刻至t = (n +1)Δt时刻曲柄角加速度an +1为:an+1 = wn+1Δt (14)以式(13)和式(14)求得曲柄转速和加速度离散点,利用MATLAB绘制成连续曲线,曲柄转速-时间曲线如图6所示,曲柄角加速度-时间曲线如图7所示。图6 曲柄转速-时间曲线Fig． 6 Curve of crank speed and time图7 曲柄角加速度-时间曲线Fig． 7 Curve of crank accleration and time4 STM32控制程序编写t = nΔt时刻的曲柄角度θ1(nΔt), t = (n + 1)Δt时刻的曲柄角度θ1((n + 1)Δt),其中, n = 0,1,2,58第12期洪 瑞等:基于牛顿迭代算法的伺服压力机控制系统研究 万方数据… ,199,电机与曲柄之间的传动比I =131．75/37．13 =3． 55,所以电机在t = nΔt时刻至t = (n +1)Δt时刻内的旋转角度为:∂n = I × Δθ1n = I × (θ1((n + 1)Δt) - θ1(nΔt))(15)式中: Δθ1 n为曲柄角度在t = nΔt时刻至t = (n +1)Δt时刻的差值。伺服驱动器预先设定一个电子齿轮比,在未改动电子齿轮比的情况下,伺服驱动器每接收10000个脉冲对应电机旋转一周,即STM32每输出10000个脉冲,电机旋转一周,所以t = nΔt时刻到t =(n +1)Δt时刻内脉冲个数numn为:numn = ■n360 × 10000 (16)t = nΔt时刻至t = (n + 1)Δt时刻内脉冲频率fren为:fren = numnΔt =■n360 × 100005 × 10-2 =■n5 × 36 × 105 =I × Δθ1n5 × 36 × 105 = 7136 × 103 × Δθ1n (17)根据公式(16)和公式(17)计算脉冲频率和个数,利用定时器TIM2实现。脉冲频率通过配置TIM2计数周期实现,关键部分代码如下:TIM_TimeBaseStructure． TIM_Period = tim2_peri-od -1; / /设置计数周期,即脉冲频率TIM_TimeBaseStructure． TIM_Prescaler = (72 -1); / /分频系数, TIM2计数频率为1 MHzTIM _ TimeBaseStructure． TIM _ CounterMode =TIM_CounterMode_Up; / /向上计数TIM_OCInitStructure． TIM_OCMode = TIM_OC-Mode_PWM1; / / PWM1输出模式TIM_OCInitStructure． TIM_Pulse = tim2_period/2 -1;脉冲个数的控制主要通过TIM2溢出中断控制,TIM2每中断溢出一次,便产生一个脉冲, pluse_number为要实现的脉冲个数, tim2_counter为每次溢出中断的计数指针,当tim2 _counter增加到和pluse_number (目标脉冲个数)相等时,意味着已实现输出指定个数脉冲,从而关闭定时器。void TIM2_IRQHandler (void) {if (TIM_GetITStatus (TIM2, TIM_IT_Update)! = RESET) {TIM_Cmd (TIM2, DISABLE); / ∗进中断之后首先关闭定时器∗ /TIM_ClearITPendingBit (TIM2, TIM_IT_Update); / ∗清除中断标志位∗ /tim2_counter + + ;if ( tim2_counter < pluse_number) {TIM_Cmd (TIM2, ENABLE);} else {tim2_counter =0;TIM_Cmd (TIM2, DISABLE);} } }电机控制主程序流程如图8所示。图8 控制程序流程图Fig． 8 Flow chart of control program最终滑块实际工作位移曲线和目标工艺曲线(拉深工艺曲线)对比如图9所示,其中实线为目标工艺曲线,虚线为滑块实际工作位移曲线,可以看出,两条曲线基本吻合。图9 滑块工作位移-时间曲线对比图Fig． 9 Comparison of curves of slider workingdisplacement and time68锻 压 技 术 第41卷万方数据5 结语提出了基于牛顿迭代法的伺服压力机非线性负载求解方案,搭建了控制系统硬件平台,利用几何方法建立四连杆机构运动学方程,利用MATLAB软件中的Simulink建立模型,通过仿真得到滑块位移-曲柄角度数据,拟合出的滑块位移与曲柄角度之间的数学方程式。基于牛顿迭代算法,在MATLAB中进行编程实现了曲柄转角的多点一次性精确求解的算法,并进一步根据求解的结果拟合出了曲柄转角-时间的函数曲线。结合伺服驱动器设置,推导出控制周期内对应的脉冲频率与个数计算方法,编写了STM32控制程序。将目标工艺曲线(拉深工艺曲线)与实际滑块工作曲线进行对比,其拟合程度较高,证明了本文控制方案的可行性。参考文献:[1] Osakada K, Mori K, Altan T, et al. Mechanical servo press tech-nology for metal forming [J]. CIRP Annals - Manufacturing Tech-nology, 2011, 60 (2): 651 -672.[2] 李麟,叶春生,莫健华.基于ARM的伺服压力机网络控制系统设计[J].锻压技术, 2013, 38 (3): 102 -106.Li L, Ye C S, Mo J H. Network control system design of servopress based on ARM [ J]. Forging & Stamping Technology,2013, 38 (3): 102 -106.[3] 周洋,叶春生,莫健华.肘杆式伺服压力机运动学仿真与工艺轨迹规划[J].锻压装备与制造技术, 2012, 47 (3): 37 -41.Zhou Y, Ye C S, Mo J H. Toggle servo press kinematics simula-tion and trajectory planning process [J]. China Metalforming E-quipment & Manufacturing Technology, 2012, 47 (3): 37 -41.[4] 雷培,刘云霞.基于MATLAB的四连杆机构运动分析[J].机械工程与自动化, 2009, (2): 74 -75.Lei P, Liu Y X. Kinematical analysis of a four-bar mechanism u-sing MATLAB [ J]. Mechanical Engineering and Automation,2009, (2): 74 -75.[5] 俞伟洋,叶春生.基于MATLAB的四连杆伺服压力机控制求解[J].锻压技术, 2015, 40 (9): 69 -73.Yu W Y, Ye C S. Research on the control solving of four-bar linkservo press based on MATLAB [J]. Forging & Stamping Technol-ogy, 2015, 40 (9): 69 -73.[6] 唐家德. MATLAB在非线性曲线拟合中的应用研究[J].智能计算机与应用, 2008, (1): 57 -59.Tang J D. Research on the application of MATLAB in nonlinearcurve fitting [J]. Intelligent Computer and Application, 2008,(1): 57 -59.[7] Dantchev S, Alhatib F. Nonlinear curve fitting for bioelectrical im-pedance data analysis: a minimum ellipsoid volume method [J].Physiological Measurement, 1999, 20 (1): N1 -9.[8] Niu Q. Numerical continuation method for nonlinear curve fittingby use of least square method [J]. Journal of Hehai University,2003, 31 (5): 597 -600.[9] 柳辉.解非线性方程的牛顿迭代法及其应用[J].重庆工学院学报:自然科学版, 2007, 21 (8): 95 -98.Liu H. Newton iterative method for solving nonlinear equations andits application [J]. Journal of Chongqing Institute of Technology:Natural Science Edition, 2007, 21 (8): 95 -98.[10] Chun C. Construction of Newton-like iteration methods for solvingnonlinear equations [ J]. Numerische Mathematik, 2006, 104(3): 297 -315.78第12期洪 瑞等:基于牛顿迭代算法的伺服压力机控制系统研究 万方数据