哈密顿系统的数学建模与动力学分析

1 引言Hamilton动力系统理论有着悠久而丰富的历史，它本身是Lagrange力学的升华与推广，从数学角度看又是一门内容精深的相空间几何学，如辛几何、辛拓扑等都源于此.近几十年来，随着纯数学理论的不断发展与计算机的普遍应用，Hamilton动力系统理论又成为当今非线性科学中极其活跃而富有魅力的研究领域.由于这类系统广泛存在于数理科学、生命科学以及社会科学的各个领域，特别是天体力学、等离子物理、航天科学以及生物工程中的很多模型都以Hamilton系统的形式出现，因此该领域的研究多年来长盛不衰.本文利用Hamilton原理推导出了Hamilton系统的正则方程.最后利用Hamilton正则方程给出一个具体物理实例的数学模型并对其进行动态模拟仿真.2 预备知识2.1 状态空间的基本概念1）状态任何一个系统在特定时刻都有一个特定的状态，系统在0t 时刻的状态是0t 时刻的一种信息量，它与此后的输入一起惟一地确定系统在0t t ≥时的行为.2）状态变量状态变量是一个完全表征系统时间域行为的的最小内部变量组. 3）状态向量设系统有n 个状态变量，用()()()12,,,n x t x t x t 表示，而且把这些状态变量看做向量()x t 的分量，则向量()x t 称为状态向量，记为()()()()12,,,Tn x t x t x t x t =⎡⎤⎣⎦.4）状态空间以状态变量()()()12,,,n x t x t x t 为轴的n 维实向量空间称为状态空间.5）状态方程描述系统状态变量与输入变量之间关系的一阶微分方程组（连续时间系统）或一阶差分方程组（离散时间系统）称为系统的状态方程，它表征了输入对内部状态的变换过程，其一般形式为：()()(),,x t f x t u t t =⎡⎤⎣⎦其中，t 是时间变量，()u t 是输入变量.6）输出方程描述系统输出量与系统状态变量和输入变量之间函数关系的代数方程称为输出方程，它表征了系统内部状态变化和输入所引起的系统输出变换，是一个变化过程.输出方程的一般形式为：()()(),,y t g x t u t t =⎡⎤⎣⎦.7）状态空间表达式状态方程与输出方程的组合称为状态空间表达式，也称动态方程，它表征一个系统完整的动态过程，其一般形式为：()()()()()(),,,,x t f x t u t t y t g x t u t t ⎧=⎡⎤⎪⎣⎦⎨=⎡⎤⎪⎣⎦⎩ 通常，对于线性定常系统，状态方程为x Ax Buy Cx Du=+⎧⎨=+⎩ 其中，()12,,Tn x x x x =表示n 维状态向量，()n n ij n nA a R ⨯⨯=∈表示系统内部状态的系数矩阵，称为系统矩阵n n A ⨯，()n r ij n rB b R ⨯⨯=∈表示输入对状态作用的矩阵，称为输入（或控制）矩阵n r B ⨯，()m n ij m nC c R ⨯⨯=∈表示输出对状态关系的矩阵，称为输出矩阵m n C ⨯，()m r ij m rD d R ⨯⨯=∈表示输入直接对输出作用的矩阵，称为直接转移矩阵m r D ⨯，也称前馈系数矩阵.A 由系统内部结构及其参数决定，体现了系统内部的特性，而B 则主要体现了系统输入的施加情况，通常情况下0D = .2.2线性定常连续系统的能控性定义2.1 设(),,n p n n x Ax Bu x R u R A R ⨯=+∈∈∈，若存在一分段连续控制向量()u t ，能在0,f t t ⎡⎤⎣⎦内，将系统从任意的初态()0x t 转移至任意终态()f x t ，则系统完全能控.定理2.1 系统完全能控的充要条件：rankSc n =其中，1,,,n Sc B AB A B -⎡⎤=⎣⎦，称为能控矩阵.2.3线性状态反馈控制律线性状态反馈控制律为u V Kx =-式中，V 是参考输入，p n K R ⨯∈称为状态反馈增益矩阵.系统动态方程变为：()()()()K K x Ax B V Kx A BK x BV A BVy Cx D V Kx C DK x DV C DV=+-=-+=+⎧⎪⎨=+-=-+=+⎪⎩ 式中，K A A BK =-，K C C DK =-，当0D =时，状态反馈系统闭环传递函数()K W s 为()()1K W s C sI A BV B -=--⎡⎤⎣⎦式中，A BV -为闭环系统的系统矩阵.以上我们简要介绍了控制系统的有关问题，现在针对单输入定常线性系统，设计其某种形式的线性定常控制律，使得闭环系统具有指定的希望的一组极点，即极点配置.2.4 极点配置考虑下述单输入线性定常系统⎩⎨⎧=+=Cx y bu Ax x （2.4.1）其中A 为n n ⨯常阵，b 和C 分别为1⨯n 和n ⨯1常阵.选取线性定常反馈控制律kx u -=，使得（2.4.1）在该控制律下的闭环系统具有指定的极点集.问题SPA[状态反馈极点配置问题] 给定矩阵n n R A ⨯∈，1⨯∈n R b 及一组共轭封闭复数i s ， i=1，2，…，n （不必互异），求取矩阵n r R K ⨯∈ 使得()i i s bK A =-λ n i ,,2,1 =对问题SPA 先考虑其解的存在性有：定义2.2 如果对于任何给定的一组共轭封闭复数i s ，n i ,,2,1 =,前述问题SPA 均有解，则称线性定常系统（2.4.1）可用状态反馈任意配置极点.下述定理给出了线性定常系统（2.4.1）利用状态反馈任意配置极点的条件. 定理2.2 定常线性系统（2.4.1）可用状态反馈任意配置极点的充要条件是系统（2.4.1）完全能控问题 对单输入系统，给定能控矩阵对[]b A 和一组期望的闭环特征值{}**2*1,,,n λλλ ，要确定n ⨯1的反馈增益矩阵k ，使成立()*i i BK A λλ=-，n i ,,2,1 =. 对于上述问题，我们有下述算法：算法2.1 [单输入系统的极点配置设计] 第一步：计算A 的特征多项式，即()0111det a s a s a s A sI n n n ++++=---第二步：计算由{}**2*1,,,n λλλ 所决定的多项式，即()()*0*11*1**1*)(a a sa s s s s a n n n n ++++=--=-- λλ 第三步：计算[]*11*11*0-----=n n a a a a a a K第四步：计算变换阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=---11][1111n n n a a a b Ab b A P 第五步：求1-=p Q第六步：所求的增益阵Q K K = .2.5 分析力学中相关的知识1） 广义坐标能够完全确定质点系位形的独立参变量，用符号 ,,21q q 表示.广义坐标是彼此独立的.其选择有一定的随意性，只需根据质点系的特点，选择那些能够惟一地确定该系统位形的参变量即可.2）广义速率在质点系中引入广义坐标之后，质点系的运动可以用广义坐标随时间的变化规律来描述，即广义速率：()k dtdq q,,2,1 ==ααα3）广义坐标变分假设在给定的运动初始条件下，某质点系的运动微分方程组的解已经求得，它的广义坐标运动方程为()t q q αα=，广义速率dtdq qαα= 于是广义坐标的全微分为 dt qdq αα = ()k ,,2,1 =α同样，广义坐标也有它的可能运动方程()t q q **αα= ()k ,,2,1 =α 比较统一瞬时广义坐标的真实运动和与其相邻的可能运动，并限定二者的差值为无限小量，即αααδq q q -=*()k ,,2,1 =α αδq 就称为广义坐标变分.4）质点系的自由度该系统独立坐标变分的数目.对完整系统它的自由度等于它的广义坐标的数目. 5）广义动量质点系的动能T 对广义速率αq的偏导数，即 ααqT p ∂∂=其中动能T 是广义坐标αq 和广义速率αq的函数. 6) 勒让德变换勒让德变换是把以n x x x ,,21为变量的函数()n x x x f ,,,21 变换成以n y y y ,,,21 为新变量的函数()n 21y …，y ,y ~f 的一种特殊变换，f ~称为f 的勒让德变换. 设有一个二次可微的函数()n x x x f ,,,21 ，且在雅可比行列式不为零，即0221212212≠∂∂∂∂∂∂∂∂∂∂nnn x f x x f x x f x f的区域内存在以下变量变换ss x fy ∂∂=()n s ,,2,1 = 定义f 的勒让德变换为()∑=-=ns s s f y x f 1n 21y …，y ,y ~于是有s sx y f=∂∂~下面给出对部分变量进行变换的情况s s x F y ∂∂=， ss y F x ∂∂=~对保留变量有rr x Fx F ~~~∂∂-=∂∂. 定理2.3 哈密顿原理从动力学普遍方程出发可推导出哈密顿原理的一般形式，即()010=+⎰dt t t W T δδ其中T δ是系统动能的变分，W δ是作用于系统的所有主动力的虚功.当作用在系统上主动力为有势时，V W δδ-=.引入哈密顿作用量⎰=tt Ldt I 10其中L 为拉格朗日函数，是系统动能与势能之差，即V T L -=.于是，对完整系统哈密顿原理可以写成常见的变分形式⎰==t t Ldt I 10δδ.3 哈密顿系统的动力学表述——哈密顿正则方程3.1 保守系统的情形拉格朗日方程是用一组关于k 个广义坐标j q 的二阶常微分方程组来描述系统的运动.方程的建立完全依赖于以()t qq j j ,, 为变量的拉格朗日函数L ，即),,,,,,,(2121t q q qq q q L L k k =.哈密顿以广义动量j p 取代广义速度j q ，以),,(t p q j j 为变量，称为哈密顿变量或正则变量.以哈密顿函数H 代替拉格朗日函数L ，用k 2个关于广义坐标j q 和广义动量j p 为变量对称整齐的一阶常微分方程组，即称为哈密顿正则方程或正则方程，以此来描述系统的运动.因本文是针对线性系统而言，故这里只给出单自由度系统的哈密顿正则方程，下面用哈密顿原理导出单自由度系统的哈密顿正则方程.首先，利用勒让德变换把以()t qq ,, 为变量的拉格朗日函数L 变换成以),,(t p q 为新变量的哈密顿函数H .显然，新变量p 代替旧变量q参与变换，而同时保留变量q 及t . 根据对原变量进行部分替换的勒让德变换，可得哈密顿函数L qp H -= 因此，拉格朗日函数H qp L -= 代入哈密顿原理，即 ()01010=-=⎰⎰dt t t H qp t t Ldt δδ 对上式进行变分运算，得0)(1=∂∂-∂∂-+⎰dt t t q qH p p H p q q p δδδδ (3.1.1) 将上式第一项改写成如下形式，即()()q pq p dtd q dt dp qp kj j j δδδδ -==∑=1 代入式(3.1.1)，有0])()[(1001=∂∂+-∂∂-+⎰dt t t q qH p p p H q t t qp δδδ (3.1.2) 因为系统在始末位置是确定的，则有0)(0=t q δ， 0)(1=t q δ (3.1.3) 于是有0])()[(10=∂∂+-∂∂-⎰dt t t q qH p p p H q δδ . (3.1.4) 根据广义动量的定义qLp ∂∂=，由部分勒让德变换可得 pHq∂∂= (3.1.5) 因此式(3.1.2)成为0])([10=∂∂+-⎰dt t t q q Hp δ对于完整系统，由于q δ可以任意取值，因此欲使上式成立，必有qHp∂∂-= (3.1.6) 联立式(3.1.5)和式(3.1.6)，即关于变量),,(t p q 的哈密顿正则方程⎪⎪⎭⎪⎪⎬⎫∂∂-=∂∂=q H pp H q. 3.2非保守系统的情形系统除有势力以外还存在非有势力作用的情形.在哈密顿原理的一般形式()010=+⎰dt t t W T δδ (3.2.1)中，系统的主动力的虚功W δ可写成如下形式：q Q V W δδδ'+-=其中，V δ-和q Q δ'分别表示有势力和非有势力的虚功.将上式代入式(3.2.1)，得dt t t q Q L dt t t q Q V T )()(1010⎰⎰'+='+-δδδδδ将H qp L -= 代入上式，并进行变分运算，得 0)(1='+∂∂-∂∂-+⎰dt t t q Q q qH p p H p q q p δδδδδ 利用式(3.1.2)和式(3.1.3)有0])()[(10='-∂∂+-∂∂-⎰dt t t q Q qH p p p H q δδ 采用与前面同样的作法，即可得到存在非有势力作用时的哈密顿正则方程pHq∂∂= Q qHp'+∂∂-= (3.2.2) 式中Q '为系统的非有势力对应于q 的广义力.4 利用哈密顿正则方程建立具体物理系统的数学模型—水平弹簧质量振动系统图 4.1 弹簧质量振动系统4.1水平弹簧质量系统的问题描述假设系统满足条件： 1) 振动无阻尼.2) 系统只能在水平方向即x 方向运动. 3) 外力()t u ，以x 的同方向为正. 要求：1) 建立弹簧质量系统的运动微分方程. 2) 求出反馈增益阵. 3) 弹簧质量系统仿真模拟. 4) 作任何有意义的讨论.4.2水平弹簧质量问题的分析解：令()t u 为输入量，y 为输出量，取弹簧等于原长0l 时，质量位置o 为x 坐标轴的原点，取y x =1为广义坐标.如势能零点取在弹簧原长位置，则系统的势能2121kx V =，因此系统的拉格朗日函数u x121212121ux kx x m V T L +-=-= . 求得广义动量1x m xLp x =∂∂=因此mp xx=1 . 计算哈密顿函数H ，并把它写成广义动量和广义坐标的函数12121212111212)2121(ux kx m p ux kx x m x p L x p H x x x -+=+--=-=求得H 后，按式(3.2.2)写出系统的正则方程mp p Hxx x =∂∂=1 u kx xHpx +-=∂∂-=1 由上二式消去x p ，得到系统运动微分方程u kx xm =+11 . 4.3 建立弹簧质量系统的数学模型令x p x =2 则有u x x u x x k x xm⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡1001001010002121121 输出方程为1x y =则弹簧质量系统的状态空间表达式⎩⎨⎧=+=Cx y bu Ax x其中 ⎥⎥⎦⎤⎢⎢⎣⎡-=o k m A 10⎥⎦⎤⎢⎣⎡=10b []01=C .5 系统闭环状态反馈控制器设计5.1系统状态反馈控制根据线性系统状态反馈控制律，设状态反馈下受控系统的输入为Kx u -= （21⨯∈R K 为反馈增益矩阵，[]21k k K =），将上式代入弹簧质量振动系统的状态空间表达式，得到弹簧质量状态反馈闭环系统的状态空间表达式⎩⎨⎧=-=Cxy x bK A x)( 其中⎥⎥⎦⎤⎢⎢⎣⎡---=-2110k kk m bk A .6 求解状态反馈增益阵由定理2.1 []⎥⎥⎦⎤⎢⎢⎣⎡==0110m Ab b s c 显然系统完全能控，故满足闭环极点可任意配置条件.取1=m ；1=k给定一组期望的闭环特征值j +-=1*1λ， j --=1*2λ1）现计算系统的特征多项式()111det det 2+=⎥⎦⎤⎢⎣⎡-=-s s s A sI 再由指定闭环极点可得希望的闭环特征多项式为()()()()2211221**++=++-+=-=∏=s s j s j s s s a i i λ 于是可求得**0011[][12]k a a a a =--=再来计算变换阵[]⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡=1001100110011011a b Ab P 并求出其逆⎥⎦⎤⎢⎣⎡==-10011P Q 从而所要确定的反馈增益阵k 即为10[12][12]01k kQ ⎡⎤===⎢⎥⎣⎦. 2）调用Matlab 函数算出的反馈增益阵见[附录1]7 动态系统的simulink仿真7.1创建Simulink系统模型首先根据弹簧质量状态反馈闭环系统的状态空间表达式，选择合适的Simulink系统模块，然后建立此系统的Simulink模型.系统的Simulink模型图见[附录3].7.2动态系统的Simulink仿真A B C D矩阵组后，,,,A B C D矩阵组表示，当输入(),,,在MATLAB中，系统状态空间用()用函数(),,,ss A B C D直接可以得到状态空间模型。