引 言
非完整约束动力系统是指约束方程中包含位移坐标对时间的导数，并且方程不可积分为有限形式的力学系统．非完整约束系统与完整约束系统的区分是Hertz等[1]开始提出的．其后，著名力学家Appell、Hamel对非完整约束动力学系统的研究做出了基础性的贡献[2-3]．在历史上，非完整约束系统的发展与几个经典例子的研究是分不开的，这些例子包括冰橇问题、滚盘问题以及滚球问题．在现代，工程技术中出现了更多的非完整系统的例子，使得非完整系统的研究又有了进一步的发展，这些发展一方面是为了解决工程技术中的某些应用；另一方面，在非完整问题的基础研究上也取得了更进一步的成果，特别是近代控制技术的发展，使得非完整约束系统的研究具备了更加广泛的意义[4]．
近四十多年来，保结构算法的研究逐步发展成为一个非常活跃的领域．在其所涵盖的众多方法之中，变分积分法是一种十分有效的方法．变分积分的思想不仅适用于Lagrange保守系统，还可以扩展到受外力和耗散影响的系统及包含约束的系统中．变分积分法基于变分原理，不同的变分原理可以构造不同的算法[5-7]．在非完整约束系统中，基于Lagrange-d’Alembert原理的一种离散形式，Cortés等先是给出了一种具有保辛格式的数值积分算法，该算法对具有对称结构的非完整约束系统，可以导出动量方程的离散形式[8];随后，在此算法的基础上，又给出了一种保能量的积分形式[9]．León等[10-11]以经典的生成函数理论为基础，针对包含线性约束的系统给出了一种可以保持系统非完整约束的积分方法，并将其应用到最优控制问题中．近年来，Mclachlan等[12]和Kobilarov等[13]也分别给出了不同形式的变分积分法，与传统算法相比，这些算法在保持系统固有的几何特性方面具有优秀的表现．上述算法都是在Lagrange体系中构造的，对于无法表示成Lagrange形式的动力系统，Gao等[14]引入了对偶向量，并基于最小作用量原理及生成函数理论，通过选择积分区间两端位移或动量为独立变量，得到了4种不同类型的求解无约束系统的保辛算法．
本文将Galerkin变分积分的思想推广到Hamilton形式下的非完整约束动力系统中．数值积分方法的构造基于对偶变量表示的Lagrange-d’Alembert原理的一种离散化形式．对于非完整约束Hamilton动力系统，利用这种离散化形式得到的数值积分方法保持了连续系统固有的对称性．保结构算法具有高阶收敛的数值特性．同时，算法能够严格满足非完整约束．
1 非完整约束系统及Lagrange-d’Alembert原理
考虑一非完整约束系统，其位置可由广义坐标q=(q1,q2,…,qd)T表示，则系统的Lagrange方程为
(1)
其中
称为Lagrange函数，是系统动能与势能之差；F称为广义力，与系统的非保守力有关．若系统受到c<d个非完整约束的作用，
(2)
则广义力F是矩阵G(q)的行向量的线性组合：
F=GT(q)λ,
(3)
其中G(q)是c×d维矩阵，λ=(λ1,λ2,…,λc)T是Lagrange乘子组成的向量．此时系统的运动方程可以表示为
(4)
方程(4)称为非完整约束系统的Lagrange-d’Alembert方程．Lagrange-d’Alembert原理指出，Lagrange-d’Alembert方程(4)由下列变分决定:
(5)
方程(5)是用Lagrange变量q和
来描述的Lagrange-d’Alembert原理．
如果引入动量
并定义Hamilton函数:
(6)
即用对偶变量q和p代替Lagrange变量q和来描述系统的运动状态，则Lagrange-d’Alembert原理可以写成
(7)
完成上式变分可以得到Hamilton正则方程：
(8)
同时，由Hamilton函数及对偶变量q和p表示的非完整约束方程可以写成
(9)
变分原理(5)和(7)均导出了非完整约束系统的运动方程，然而变分原理(5)由一类变量表示，变分原理(7)由对偶变量表示，二者在构造数值算法时是不同的．应看到，两种形式的Lagrange-d’Alembert原理都不是标准的变分原理，非完整约束方程都不能由变分原理导出．然而，大量文献表明[8-13]，对于方程(5)，即使其不是标准的变分原理，依然可以结合变分积分的思想来构造数值积分方法．基于这一点，本文将利用对偶变量表示的变分原理(7)，研究求解非完整约束系统的保结构算法．
2 非完整约束系统保结构算法
回顾无约束系统中Galerkin变分积分法的构造过程[5,15]，算法的构造基于对作用量的离散，即选取一个有限维的多项式空间来近似广义位移，同时选取适当的数值积分方法来近似作用量中的积分．拓展到非完整约束系统中，由于Lagrange-d’Alembert原理并非标准的变分原理，没有与之对应的作用量，本文考虑直接对变分方程进行离散，从而得到数值算法．
为此，在区间[0,η]上，首先选取适当的多项式，对Lagrange-d’Alembert原理(7)中所包含的广义位移q、广义动量p以及Lagrange乘子λ进行近似．以广义位移为例，选取s+1个Lobatto点0=d0<d1<…<ds-1<ds=1及控制点q0=q(0),q1,q2,…,qs-1,qs=q(η)，从而唯一地确定一个s阶多项式q(t)，使得对于所有的Lobatto点和控制点满足q(dvη)=qv，v=0,1,2,…,s，其中qv为广义位移在dvη时刻的值．若用lv(t)表示Lagrange多项式，则广义位移的近似多项式q(t)可以表示为
(10)
类似地，可以给出广义动量和Lagrange乘子的近似多项式：
(11)
其中p0=p(0),p1,p2,…,ps-1,ps=p(η)和λ0=λ(0),λ1,λ2,…,λs-1,λs=λ(η)分别为广义动量和Lagrange乘子的控制点．
然后，将式(10)和(11)表示的广义坐标、广义动量以及Lagrange乘子的近似多项式代入变分原理(7)中，同时选取s+1个点的Lobatto积分法近似变分原理中的积分项，记τi和wi分别为Lobatto积分法在[-1,1]区间上的积分点及权重，i=0,1,…,s,则变分原理(7)离散为
(12)
其中
(13)
考虑以时间区间两端位移为独立变量的生成函数Sd(q(0),q(η))，其对应的正则变换如下：
(14)
完成方程(12)的变分并结合正则变换(14)，可得到如下方程：
(15)
其中
将区间内的非独立变量
和
写成独立变量q0和qs函数的形式．
注意到q0，p0和λ0为已知量，方程组(15)中包含2(s+1)d+sc个未知数，但仅包含2(s+1)d个方程，因此算法还需要sc个方程．本文令非完整约束在区间端点处以及区间内部的控制点处均精确满足，从而得到sc个方程:
(16)
求解非线性代数方程组(15)和(16)便可得到非完整约束系统的解，本文选取Newton迭代法求解．至此，基于Galerkin变分积分法和对偶变量p和q表示的Lagrange-d’Alembert原理(7)，本文得到了一种求解非完整约束Hamilton动力系统的数值积分方法．
3 对称性证明
非完整约束Hamilton动力系统为可逆系统，对称算法在求解可逆系统时，具有明显的数值优势[16]．对于一个数值方法y1=ψη(y0)，若其满足
(17)
即算法在交换y0↔y1和η↔-η后保持不变，称该算法为对称的或时间可逆的[16]．因此，对于本文所构造的算法，如果用-η替换η，用qs替换q0，ps替换
替换
替换
替换
而方程组(15)和(16)保持不变，则说明算法是对称的．下文将证明，在实施上述替换之后，方程Fv将转化为Fs-v，Ev将转化为Es-v，Dv将转化为Ds-v，方程组整体保持不变，从而证明了算法的对称性．
将式(13)代入F0和Fs中得到展开式:
(18)
(19)
其中
(20)
由于lj(t)为Lagrange多项式，τi和wi分别为Lobatto积分法在[-1，1]区间上的积分点及权重，容易得到如下关系：
lj,i=ls-j,s-i, l′j,i=-l′s-j,s-i, τi=-τs-i, wi=ws-i．
(21)
将式(21)代入式(18)，同时用-η替换η，用qs替换q0，ps替换
替换
替换
替换
得到
(22)
容易验证式(22)等价于式(19)，即方程F0转化为方程Fs．类似地，我们可以证明方程Fv将转化为方程Fs-v，v=1,2,…,s-1，Ev将转化为Es-v，v=0,1,…,s．
注意到，初始值q0和p0为已知量，自动满足约束条件，若记初始点处的约束方程为
(23)
对于方程(23)和(16)，如果用qs替换q0，ps替换替换替换
很容易验证方程Dv将转化为方程Ds-v．综上可证，本文算法具有对称性．
4 数 值 算 例
算例1 首先考虑非完整约束系统的经典算例——冰橇问题[4]．设一冰刀在倾斜平面上运动，系统的Lagrange函数为
(24)
非完整约束方程为
(25)
该系统可求得解析解，如果初始条件为
(26)
则解析解为
(27)
系统参数分别取m=1，J=1，g=10，α=π/12，ω=1，采用本文算法对该算例进行仿真．该算例具有解析解，可以给出如下的状态向量相对误差ev及全局误差e的定义：
(28)
其中va为状态向量的解析解，vk为kη时刻算法得到的数值解．图1给出了使用不同插值阶数s=1,2,…,8和积分步长η=1/2n，n=1,2,…,7时，本文算法所得到的状态向量全局误差，曲线斜率表示算法的收敛阶数．从图1可以看出，当插值阶数s给定后，算法的收敛阶数O满足
(29)
为观察算法长时间仿真性质，图2给出了仿真总时长为2 000，插值阶数s=3，积分步长η=1/8时的仿真结果．图2(a)为系统状态向量的相对误差曲线，图2(b)为Hamilton函数相对误差曲线．对于保守系统，Hamilton函数是守恒量，其相对误差为eH=(Hk-H0)/H0，其中H0为Hamilton函数的初始值，Hk为kη时刻算法得到的数值解．图2(b)表明，Hamilton函数相对误差始终在明确的幅值范围内波动，说明算法具有良好的能量属性．图2(c)为非完整约束的仿真值，数值结果达到机器精度，说明本文算法严格满足非完整约束．
图1 插值阶数s和积分步长η取不同值时，状态向量全局误差
Fig. 1 The global errors of state vectors with different orders of interpolation s and step sizes η
算例2 本算例研究蛇板的运动，其属于一类典型的带有非完整约束的多体系统．在Ostrowski等[17]建立的蛇板系统的动力学模型中，系统的Lagrange函数为
(30)
非完整约束方程为
(31)
(a) 状态向量相对误差 (b) Hamilton函数相对误差 (c) 非完整约束
(a) The relative error of the state vector (b) The relative error of the Hamiltonian function (c) The nonholonomic the state vector constraint
图2 冰刀长时间仿真结果
Fig. 2 Numerical results of the skate for long-time simulation
本文中，系统参数为m=6,l=0.2,J=0.016,Jr=0.072,Jw=0.001 3．分别选取如下3种初值条件进行研究．第一种初值条件：
(32)
第二种初值条件：
(33)
第三种初值条件：
(34)
(a) 第一种初值条件(b) 第二种初值条件 (c) 第三种初值条件
(a) The 1st initial condition (b) The 2nd initial condition (c) The 3rd initial condition
图3 蛇板运动轨迹
Fig. 3 The orbits of the snakeboard
图3分别给出了3种初值条件下蛇板的运动轨迹．采用插值阶数s=4，积分步长η=1/16的算法进行仿真，积分时长2 000．图3(a)表明，仿真结果得到的位移轨迹与实际运动情况相吻合．且经过长时间仿真后，蛇板仍保持匀速直线运动，说明算法具有良好的稳定性．在第二种初值条件下，α取不同值会产生不同的运动轨迹．针对这一初值条件，采用两种不同参数的算法进行计算，对比结果的差异．图3(b)给出了α取不同值时两种参数算法求得的运动轨迹，图中实线表示s=4,η=1/16的算法求得的结果，“∘”表示s=1，η=1/2的算法求得的结果．可以看出，两种参数的算法求得的位移轨迹完全一致，说明本文算法即使采用较低的收敛阶数和较大的积分步长，仿真结果依然收敛．图4为第二种初值条件下，s=1,η=1/2时，算法求得的非完整约束．图4表明，即使采用最低的插值阶数和较大的积分步长，算法依然严格满足非完整约束．最后，图5给出了第三种初值条件下，s=4，η=1/16，积分时长为2 000的仿真结果．从图中可以看出，在复杂的运动状态下，算法依然保持良好的长时间积分性质．
(a) α=π/8 (b) α=π/16 (c) α=π/32
图4 s=1，η=1/2时，蛇板的非完整约束
Fig. 4 The nonholonomic constraints of the snakeboard for s=1,η=1/2
图5 第三种初值条件下，蛇板长时间仿真结果
Fig. 5 Numerical results of the snakeboard for long-time simulation with the 3rd initial condition
5 结 论
本文利用Galerkin变分积分的思想，以Hamilton函数和对偶变量表示的Lagrange-d’Alembert原理为基础，给出了一类求解非完整约束Hamilton动力系统的保结构算法．算法具有高阶收敛性、对称性，并且能够严格满足非完整约束．算法构造时，若使用s阶Lagrange多项式对广义位移、广义动量及Lagrange乘子进行近似，使用s+1个点的Lobatto积分法对变分原理中的积分进行近似，那么当s为奇数时，算法的收敛阶数为O=s+1;当s为偶数时，算法的收敛阶数为O=s．
[1] HERTZ A, GARBASSO A. Die prinzipien der mechanik in neuem zusammenhang dargestellt[J]. Il Nuovo Cimento (1895—1900), 1895, 1(1): 40-59.
[2] APPELL P. Traité de Mécanique Rationnelle[M]. Gauthier-Villars, 1924.
[3] HAMEL G. Theoretische Mechanik[M]. Berlin: Springer, 1978.
[4] BLOCH A M. Nonholonomic Mechanics and Control[M]. Berlin: Springer, 2015.
[5] MARSDEN J E, WEST M. Discrete mechanics and variational integrators[J]. Acta Numerica, 2001, 10(1): 357-514.
[6] LEW A, MARSDEN J E, ORTIZ M, et al. Variational time integrators[J]. International Journal for Numerical Methods in Engineering, 2010, 60(1): 153-212.
[7] HAIRER E, WANNER G, LUBICH C. Geometric Numerical Integration[M]. Berlin: Springer, 2002.
[8] CORTÉS J, MART
NEZ S. Non-holonomic integrators[J]. Nonlinearity, 2001, 14(5): 1365-1392.
[9] CORTÉS J. Energy conserving nonholonomic integrators[J]. Discrete & Continuous Dynamical Systems, 2003, 2003(S): 189-199.
[10] DE LE
N M, DE DIEGO D M, SANTAMARIA-MERINO A. Geometric integrators and nonholonomic mechanics[J]. Journal of Mathematical Physics, 2002, 45(3): 1042.
[11] DE LEN M, DE DIEGO D M, SANTAMARIA-MERINO A. Geometric numerical integration of nonholonomic systems and optimal control problems[J]. European Journal of Control, 2004, 10(5): 515-521.
[12] MCLACHLAN R, PERLMUTTER M. Integrators for nonholonomic mechanical systems[J]. Journal of Nonlinear Science, 2006, 16(4): 283-328.
[13] KOBILAROV M, MARSDEN J E, SUKHATME G S. Geometric discretization of nonholonomic systems with symmetries[J]. Discrete and Continuous Dynamical Systems(Series S), 2010, 3(1): 61-84.
[14] GAO Q, TAN S J, ZHANG H W, et al. Symplectic algorithms based on the principle of least action and generating functions[J]. International Journal for Numerical Methods in Engineering, 2012, 89(4): 438-508.
[15] HALL J, LEOK M. Spectral variational integrators[J]. Numerische Mathematik, 2012, 130(4): 681-740.
[16] HAIRER E, LUBICH C, WANNER G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[M]. Berlin: Springer, 2006.
[17] OSTROWSKI J, LEWIS A, MURRAY R, et al. Nonholonomic mechanics and locomotion: the snakeboard example[C]//Paper Presented at the Proceedings of the 1994 IEEE International Conference on Robotics and Automation. Piscataway, USA, 1994.