最小二乘法在高斯拘束函數(shù)求解中的應(yīng)用*

袁萍萍戈新生

（1.北京信息科技大學(xué)機(jī)電工程學(xué)院，北京 100192）（2.北京信息科技大學(xué)理學(xué)院，北京 100192）

最小二乘法在高斯拘束函數(shù)求解中的應(yīng)用*

袁萍萍1?戈新生2

（1.北京信息科技大學(xué)機(jī)電工程學(xué)院，北京 100192）（2.北京信息科技大學(xué)理學(xué)院，北京 100192）

高斯最小拘束原理是一種典型的微分變分原理，以加速度為變量，通過尋求拘束函數(shù)極值的變分方法直接得出系統(tǒng)的運(yùn)動(dòng)規(guī)律.目前，國內(nèi)常用的求解高斯拘束函數(shù)的方法為拉格朗日乘子法，通過引入拉格朗日乘子將高斯拘束函數(shù)的條件極值問題轉(zhuǎn)化為帶有拉格朗日乘子的無條件極值問題，這種求解方法會(huì)增加未知變量的個(gè)數(shù).為減少變量個(gè)數(shù)，進(jìn)一步提高運(yùn)算效率，文章首先對高斯拘束函數(shù)進(jìn)行簡單變形，引入加速度形式的約束方程將高斯拘束函數(shù)化為最小二乘形式，直接運(yùn)用最小二乘法導(dǎo)出使高斯拘束函數(shù)取極小值時(shí)系統(tǒng)真實(shí)加速度的表達(dá)式.最后通過對曲柄滑塊機(jī)構(gòu)正動(dòng)力學(xué)問題的分析和計(jì)算，驗(yàn)證了該方法的有效性.

多體系統(tǒng)， 動(dòng)力學(xué)建模， 高斯最小拘束原理， 最小二乘法， 曲柄滑塊機(jī)構(gòu)

引言

多體系統(tǒng)動(dòng)力學(xué)廣泛應(yīng)用于航天、機(jī)械、兵器、車輛等領(lǐng)域，為分析多體系統(tǒng)的運(yùn)動(dòng)性能及優(yōu)化設(shè)計(jì)提供了有力的理論支撐.多體系統(tǒng)動(dòng)力學(xué)的研究包括系統(tǒng)動(dòng)力學(xué)模型的建立和提高模型求解的速度、穩(wěn)定性和精度等內(nèi)容.栗鵬和姚文麗［1］以共線三球鏈的碰撞動(dòng)力學(xué)為例，建立其碰撞的Hertz模型并進(jìn)一步研究了求解的數(shù)值算法.目前，常用的解決多體系統(tǒng)動(dòng)力學(xué)問題的方法有：牛頓-歐拉方法、羅伯遜-維登伯格方法、拉格朗日方法和高斯最小拘束原理等.高斯最小拘束原理屬于微分變分原理，變分的力學(xué)原理是將系統(tǒng)的真實(shí)運(yùn)動(dòng)與可能發(fā)生的運(yùn)動(dòng)進(jìn)行區(qū)分并加以比較.高斯最小拘束函數(shù)以系統(tǒng)的加速度為變量，在滿足約束方程的系統(tǒng)所有可能加速度中，真實(shí)加速度使拘束函數(shù)取極小值［2］.運(yùn)用高斯最小拘束原理可以直接得出系統(tǒng)的運(yùn)動(dòng)規(guī)律.近年來，陳小平、馮祖仁等［3］將高斯最小拘束原理應(yīng)用于端桿加速度受約束的機(jī)器人動(dòng)力學(xué)分析和建模中，得出機(jī)器人關(guān)節(jié)加速度的線性遞推公式.史躍東和王德石［4］運(yùn)用高斯最小拘束原理建立了火炮系統(tǒng)的動(dòng)力學(xué)模型，進(jìn)行了火炮射擊過程的數(shù)值仿真和振動(dòng)分析，并確定了影響炮口振動(dòng)的主要因素以及火炮系統(tǒng)的固有振動(dòng)頻率.董龍雷、閆桂榮等［5］結(jié)合有限段法和高斯最小拘束原理對一類剛?cè)峄旌舷到y(tǒng)進(jìn)行描述，比較完整的反映了系統(tǒng)的運(yùn)動(dòng)學(xué)特征.姚文麗和戴葆青［6］從相對質(zhì)點(diǎn)形式的變分形式推導(dǎo)出廣義坐標(biāo)形式的高斯最小拘束原理，并給出了非理想約束、單邊約束及剛體碰撞情形下的高斯最小拘束函數(shù)表達(dá)式.劉延柱、薛紜［7］建立了基于高斯原理的Cosserat彈性桿的分析力學(xué)模型，導(dǎo)出了拘束函數(shù)的普遍形式，將彈性桿的相應(yīng)約束對可能運(yùn)動(dòng)施加限制，避免彈性桿空間自相侵占情況的發(fā)生.Kalaba［8］運(yùn)用拉格朗日方程的一般形式推出高斯最小拘束函數(shù)表達(dá)式，并導(dǎo)出滿足約束條件的系統(tǒng)加速度表達(dá)式.Herman和Oussama［9］將高斯最小拘束原理應(yīng)用于冗余系統(tǒng)中，導(dǎo)出冗余機(jī)器人的動(dòng)力學(xué)方程.Fan和Kalaba［10］將高斯最小拘束原理運(yùn)用到非理想約束系統(tǒng)中，分析非理想約束力，如摩擦力，對系統(tǒng)運(yùn)動(dòng)精度的影響.

目前，論文中常用的求解高斯拘束函數(shù)的拉格朗日乘子法，需引入h個(gè)拉格朗日乘子λi（i＝1，2，…，h），其中h為約束方程的個(gè)數(shù)，雖然此方法可以將高斯拘束函數(shù)的條件極值問題轉(zhuǎn)化為無條件極值問題，但是相應(yīng)的增加了h個(gè)未知變量，降低了計(jì)算效率，為避免引入多余的未知變量，文章試圖尋求一種新的求解高斯拘束函數(shù)的方法.

通過簡單的變形運(yùn)算，將求高斯拘束函數(shù)最小值問題變形為最小二乘問題.根據(jù)系統(tǒng)的真實(shí)加速度滿足約束方程且使拘束函數(shù)取極小值，運(yùn)用最小二乘法導(dǎo)出系統(tǒng)真實(shí)加速度的表達(dá)式，對系統(tǒng)的真實(shí)加速度進(jìn)行積分運(yùn)算便可得到系統(tǒng)的真實(shí)運(yùn)動(dòng)規(guī)律.最后運(yùn)用推導(dǎo)出的系統(tǒng)真實(shí)加速度的顯示表達(dá)式，對曲柄滑塊機(jī)構(gòu)進(jìn)行動(dòng)力學(xué)分析及計(jì)算，得到其運(yùn)動(dòng)規(guī)律，驗(yàn)證該求解方法的有效性.

1 高斯最小拘束函數(shù)

高斯最小拘束函數(shù)以加速度為變量，在任意時(shí)刻t，系統(tǒng)的真實(shí)運(yùn)動(dòng)與符合系統(tǒng)約束的所有可能發(fā)生的運(yùn)動(dòng)相比較，系統(tǒng)的真實(shí)加速度使其拘束函數(shù)取極小值.

高斯最小拘束函數(shù)G的表達(dá)式為［8］

式（2）為高斯最小拘束函數(shù)的另一種表達(dá)形式.式（1）和（2）考慮了系統(tǒng)所受理想約束的條件.當(dāng)Qc為非理想廣義約束力時(shí)，可將非理想約束力視為系統(tǒng)主動(dòng)力的一部分，此時(shí)高斯最小拘束函數(shù)為

2 系統(tǒng)的約束方程

運(yùn)用高斯最小拘束原理求解多體系統(tǒng)的運(yùn)動(dòng)規(guī)律，其系統(tǒng)的約束方程可以是完整的也可以是非完整的.設(shè)系統(tǒng)受到d個(gè)完整約束、g個(gè)非完整約束，其形式可寫為

將式（4）對時(shí)間求兩次微分和式（5）對時(shí)間求一次微分，所得結(jié)果可寫為矩陣形式

3 求解高斯最小拘束函數(shù)中加速度

4 仿真算例

以曲柄滑塊機(jī)構(gòu)為例，如圖（1）所示.圖中o1，o2和o3分別為曲柄、連桿和滑塊的質(zhì)心，曲柄、連桿和滑塊的連體坐標(biāo)系建立在其質(zhì)心處，慣性坐標(biāo)系建立在點(diǎn)o處，曲柄和連桿的連體坐標(biāo)系和慣性坐標(biāo)系之間的夾角分別為θ1，θ2.系統(tǒng)的廣義坐標(biāo)為：q＝［x1，y1，θ1，x2，y2，θ2，x3］T.曲柄在轉(zhuǎn)動(dòng)副o處施加有轉(zhuǎn)矩τ＝-0.5N·m.滑塊末端連接有彈簧阻尼器，剛度系數(shù)為k＝5N/m，阻尼系數(shù)為c＝1N·s/m.當(dāng)θ1＝θ2＝0時(shí)，彈簧處于平衡位置.設(shè)m1，m2和m3分別為曲柄、連桿和滑塊的質(zhì)量，m1＝0.37kg，m2＝0.77kg，m3＝0.1kg；l1和l2分別為曲柄和連桿的長度，l1＝0.15m，l2＝0.56m；g為重力加速度.

圖1 曲滑塊機(jī)構(gòu)Fig.1 Crank slidermechanism

曲柄滑塊機(jī)構(gòu)的廣義質(zhì)量矩陣M，廣義力向量Q分別為

系統(tǒng)的約束方程為

設(shè)系統(tǒng)的初始條件分別為

將式（17）、（18）、（19）、（20）代入式（3），得到曲柄滑塊機(jī)構(gòu)的拘束函數(shù)表達(dá)式.引入約束方程φ，通過式（16）求解得出曲柄滑塊機(jī)構(gòu)的位移和速度變化規(guī)律圖.

圖1～5分別為曲柄質(zhì)心在X，Y方向上的位移、速度變化曲線.圖6～7為曲柄轉(zhuǎn)角與曲柄轉(zhuǎn)動(dòng)角速度的變化曲線.圖8～9為滑塊質(zhì)心在X方向的位移、速度變化曲線.曲柄滑塊機(jī)構(gòu)中存在阻尼為耗散機(jī)構(gòu)，滑塊的運(yùn)動(dòng)幅度和速度會(huì)隨時(shí)間逐漸減小，機(jī)構(gòu)最終在一平衡狀態(tài)停止運(yùn)動(dòng).

圖2 曲柄質(zhì)心X軸方向位移Fig.2 Deformation of the crank shaft in X axis direction

圖3 曲柄質(zhì)心X軸方向速度Fig.3 Velocity of Crank center in X direction

圖4 曲柄質(zhì)心Y軸方向位移Fig.4 Displacement of crank Center in Y direction

圖6 曲柄轉(zhuǎn)角隨時(shí)間的變化規(guī)律Fig.6 Time history of the crank angle

圖7 曲柄轉(zhuǎn)動(dòng)角速度Fig.7 Angular velocity of Crank

圖8 滑塊X軸方向的位移Fig.8 Displacement of the slider in X direction

圖9 滑塊X軸方向的速度Fig.9 Velocity of the slider in X direction

5 小結(jié)

高斯最小拘束原理不直接描述系統(tǒng)的機(jī)械運(yùn)動(dòng)規(guī)律，而是把真實(shí)發(fā)生的運(yùn)動(dòng)與可能發(fā)生的運(yùn)動(dòng)加以比較，在相同條件下所發(fā)生的很多的可能運(yùn)動(dòng)中指出真實(shí)運(yùn)動(dòng)所應(yīng)滿足的極值條件，用高斯最小拘束原理求解多體系統(tǒng)動(dòng)力學(xué)問題，主要優(yōu)點(diǎn)是可以利用各種有效的數(shù)學(xué)規(guī)劃方法尋求高斯拘束函數(shù)的極值，且不論是樹形的還是非樹形的系統(tǒng)，都可以用同樣的方法處理.文章利用最小二乘法對高斯拘束函數(shù)進(jìn)行求解，直接得到了系統(tǒng)真實(shí)運(yùn)動(dòng)的加速度表達(dá)式，且不需要引入其它輔助變量.從曲柄滑塊機(jī)構(gòu)的仿真結(jié)果可以看出，系統(tǒng)的數(shù)學(xué)模型推導(dǎo)簡單，形式統(tǒng)一，而且能夠得到多體系統(tǒng)運(yùn)動(dòng)規(guī)律的精確解，計(jì)算效率高，方法簡單，易于程式化.本文提供了一種統(tǒng)一而有效的多體系統(tǒng)動(dòng)力學(xué)分析與計(jì)算方法.

1 栗鵬，姚文麗.共線三球鏈問題的碰撞動(dòng)力學(xué)研究.動(dòng)力學(xué)與控制學(xué)報(bào)，2013，11（4）：301～305（Li P，YaoW L. Study on collinear collision dynamics of the three-ball chain.Journal of Dynamicsand Control，2013，11（4）：301～305（in Chinese））

2 劉延柱.高等動(dòng)力學(xué).北京：高等教育出版社，2001（Liu Y Z.Higher dynamics.Beijing：Higher Education Press，2001（in Chinese））

3 陳小平，馮祖仁，胡保生.一種新的基于spatial高斯最小拘束原理的正動(dòng)力學(xué)算法.機(jī)器人，1992，14（2）：45～51（Chen X P，F(xiàn)eng Z R，Hu B S.A new positive dynamic algorithm based on the principle of spatial Gauss. Robot，1992，14（2）：45～51（in Chinese））

4 史躍東，王德石.高斯最小拘束原理在火炮振動(dòng)分析中的應(yīng)用.火炮發(fā)射與控制學(xué)報(bào)，2009（4）：26～30（Shi Y J，Wang D S.Gaussianminimum binding principle in Gun Vibration Analysis.Gun Launch＆Control，2009（4）：26～30（in Chinese））

5 董龍雷，閆桂榮等.高斯最小拘束原理在一類剛?cè)狁詈舷到y(tǒng)分析中的應(yīng)用.兵工學(xué)報(bào)，2001，22（3）：347～351（Dong L L，Yan GR，Du Y T.Gaussianminimum binding principle in a class of rigid coupling system analysis.Industry and Technololgy，2001，22（3）：347～351（in Chinese））

6 姚文莉，戴葆青.廣義坐標(biāo)形式的高斯最小拘束原理及其推廣.力學(xué)與實(shí)踐，2014，36（6）：779～785（Yao W L，Dai B Q.Generalized coordinates Gaussian minimum binding principle and its promotion.Mechanics in Engineering，2014，36（6）：779～785（in Chinese））

7 劉延柱，薛紜.基于高斯原理的Cosserat彈性桿動(dòng)力學(xué)模型.物理學(xué)報(bào)，2015，64（4）：1～5（Liu Y Z，Xue Y.Dynamic model of Cosserat elastic rod based on gauss principle.Journal of Physics，2015，64（4）：1～5（in Chinese））

8 Kalaba R E，Udwadia F E.Equations ofmotion for nonholonomic，constrained dynamical systems via gauss′s principle. Journal of Applied Mechanics，1993，60（3）：662～668

9 Herman B，Oussama K.Gauss′s principle and the dynamicsof redundant and constrained manipulators.International Conference on Robotics＆Automation，2000，3：2563～2568

10 Fan Y Y，Kalaba R E.Reflections on the gauss principle of least constraint.Journal of Optimization Theory and Applications，2005，127（3）：475～484

11戴中林.廣義逆矩陣A-的計(jì)算方法及應(yīng)用.大學(xué)數(shù)學(xué)，2014，30（6）：56～59（Dai Z L.The calculation method of generalized inversematrix A-and its application.University Mathematics，2014，30（6）：56～59（in Chinese））

12李慶揚(yáng)，王能超等.數(shù)值分析.北京：清華大學(xué)出版社，2001（LiQ Y，Wang N C，Yi D Y.Numerical analysis. Beijing：Tsinghua University Press，2001（in Chinese） ）

APPLICATION OF LEAST SUQARE METHOD IN SOLVING GAUSS CONSTRAINT FUNCTION*

Yuan Pingping1?Ge Xinsheng2
（1.School of Mechanical and Electrical Engineering，Beijing Information Science and Technology University，Beijing 100192，China）
（2.School of Applied Science，Beijing Information Science and Technology University，Beijing 100192，China）

Gauss′sminimum constraint principle is a typical differential variation principle，where the acceleration is the variable，and themotion law of the system can be obtained directly by the variationmethod of seeking the extreme of the constraint function.At present，Lagrange multiplier method is widely used for solving the Gauss constraint function.Through the introduction of Lagrangemultipliers，the conditional extreme value problem of Gauss constraint function is transformed to the unconditional extreme problem.However，thismethod increase the number of unknown variables.In order to reduce the number of variables，it need further study to improve the operation efficiency.In this paper，the deformation on Gauss constraint function is firstly simplified. The constraint equations of the acceleration form are also introduced into the Gauss constraint equations to get the least square form of equations.The least squaremethod is then used to export the expression of the true acceleration of the system，which canmake the Gauss binding function to take theminimum value.In the end，the validity of themethod is verified by analyzing and calculating the normal dynamics of the crank slidermechanism.

multi-body system， dynamic modeling， Gauss′sminimum constraint principle， Least square method， Crank slidermechanism

10.6052/1672-6553-2016-003

2015-10-08收到第1稿，2015-12-08收到修改稿.

*國家自然科學(xué)基金資助項(xiàng)目（11472058）

?通訊作者E-mail：pingping861206＠163.com

Received 8 October 2015，revised 8 December 2015.

*The project supported by the National Natural Science Foundation of China（11472058）

?Corresponding author E-mail：pingping861206＠163.com