Hopf分岔是非线性动力系统中相对比较简单而又重要的分岔问题，属于一种局部动态分岔^{[1, 2]}，随着分岔参数的变化，系统在非双曲平衡点处从平衡点突然分岔出极限环.对于非线性系统Hopf分岔的研究，一般需要解决如下3个方面的问题：1)Hopf分岔的存在性，即系统是否存在周期解，这是Hopf分岔研究中需要解决的基本问题；2)Hopf分岔的方向，即出现分岔的参数范围；3)Hopf分岔的稳定性，即存在周期解时，其稳定性如何(超临界或亚临界).

Hopf分岔问题的提出及其理论研究均具有深刻的应用背景^[3-6]，它通常与系统自激振荡有密切联系.如在水轮机非线性系统中，即使无任何外部激励，也可能会出现某一固定频率和振幅的自激振荡，还有可能出现频率和振幅不相同的自激振荡，甚至出现更复杂的非线性振荡以至于混沌现象^[7].然而，到目前为止，对水轮机调节系统中可能存在的分岔与振荡问题却一直研究较少，水轮机调节系统中存在着复杂分岔现象也是理论工作者和工程技术人员所关心的问题.对水轮机调节系统而言，分岔现象是一种严重的失稳现象，通常与系统的自激振荡有密切联系，因此，在控制系统设计中必须给予足够的重视^[8-9].

双Hopf是一类非常重要的余维Hopf分岔现象，此时，从平衡点分岔出来的不再是周期轨道，而可能是拟(准)周期轨道，在适当选取的三维相空间中形成轮胎面或指环形的吸引域.文献[10-13]运用多尺度方法计算系统的规范形，借助相关计算软件，对一些较典型的机械及电路系统中的双Hopf分岔进行了分析计算.但针对水轮机调节系统的双Hopf分岔分析研究却很少，因此，揭示水轮机调节系统中存在的双Hopf分岔现象和发生机理很有意义.

鉴于此，本文首先给出水轮机调节系统刚性水击和弹性水击时对应的四维和六维模型，并说明了系统Hopf存在的条件及分岔方向的计算，再根据系统特征多项式系数推导出双Hopf分岔存在的判别和计算方法，最后通过仿真和数值计算对理论分析进行验证.

1 Hopf分岔存在性判定

定理1：考虑C^k(k≥0)向量场f：U×V∈Rⁿ上的含参数的非线性动力系统：

其中，v为系统的参数，称为分岔参数或控制变量.设向量场f∈C²，(x₀，v₀)∈U×V是平衡点，若下列条件成立，则在平衡点(x₀，v₀)处出现极限环：

1) f(x₀，v₀)=0；

2) Df_x(x₀，v₀)有一对纯虚特征值λ(v₀)=±iω₀，其余n-2个特征根的实部均不为0；

3) a′(v₀)=d/dvReλ(v₀)≠0.

定理1称为Hopf分岔定理，条件2表明Hopf分岔点是非双曲平衡点，Hopf分岔现象发生在一对共轭特征值穿越虚轴时，常用于作Hopf分岔现象的判据；条件3称为横截系数，主要用于判断Hopf分岔的方向.

但定理1并未给出Hopf分岔的稳定性判定，因此需要进一步计算与Hopf分岔方向及稳定性相关的特征量.一般来说，Hopf分岔存在性问题只需要用到原非线性系统线性化部分的信息，而稳定性问题的解决则需要用到非线性部分的信息，与存在性问题相比，稳定性问题更为复杂，从本质上讲其原因在于存在性问题仅涉及到解的拓扑结构，而稳定性问题则与解的度量性质有关，将会涉及到更复杂的计算.

描述Hopf分岔方向及稳定性的特征量一般包括：1)系统Jacobi矩阵特征值实部的导数a′(v₀)，即横截系数；2)曲率系数β值.根据分岔的横截系数(或规范形系数)和曲率系数可将Hopf分岔分为超临界和亚临界两种，超临界Hopf分岔对应于稳定的振荡或系统状态相空间中的稳定极限环；而亚临界Hopf分岔对应于不稳定的振荡或系统状态相空间中的不稳定极限环，将对系统的稳定运行产生不利的影响.

横截系数a′(v₀)和曲率系数β与Hopf分岔的稳定性及分岔方向之间的关系可归结为下述定理.

定理2^[14]：当系统存在Hopf分岔，若β < 0，表明系统发生超临界的Hopf分岔，系统的周期解是稳定的；当a′(v₀)>0，稳定的极限环出现在v₀>0一侧，而在v₀ < 0的一侧系统是渐近稳定的；而当a′(v₀) < 0，极限环出现在v₀ < 0一侧，在v₀>0的一侧系统是渐近稳定的.若β>0，则表明系统的Hopf分岔是亚临界的；a′(v₀)>0时，极限环出现在v₀ < 0一侧充分小的邻域内，a′(v₀) < 0时，极限环出现在v₀>0一侧充分小的邻域内.

由定理1，当系统中存在非双曲平衡点，且其Jacobi矩阵中出现一对纯虚特征值，其余特征值均有负实部时，则系统发生Hopf分岔.而在四维及四维以上的系统中，当Jacobi矩阵出现多对纯虚特征值时，则系统发生的是多Hopf分岔，其中以双Hopf分岔(Hopf-Hopf分岔)较为常见，即两对复共轭特征值同时穿越虚轴的现象，其判别方法在后文中根据具体系统将详细介绍.如图 1(a)所示为两对复共轭特征值同时从复平面左半侧进入右半平面的情况，表明此时系统从稳定状态通过分岔失稳；图 1(b)为两对复共轭特征值中的一对从复平面左半侧进入右半平面的同时，另一对复共轭特征值从复平面右半侧进入左半平面，表明系统在经历双Hopf分岔前后均处于不稳定状态.

1.1 四维水轮机调节系统

对于一般的非线性系统，特别是维数较高时，直接应用定理1来判断Hopf分岔的存在性并不方便.实践中需要更直接的方法用于Hopf分岔存在性的判断.对于四维非线性系统，定理3给出了系统存在Hopf分岔的条件，并可同时计算出发生Hopf分岔时系统分岔值.

定理3：非线性系统y=f(x, v)中，x∈R⁴，v∈R¹为分岔参数，x=0为系统的平衡点，在平衡点时的系统Jacobi矩阵特征多项式可表示为

若当v=0时，以下各条件成立：

1) 系数p_i>0(i=0, 1, 2, 3)，且满足：

2) 其中，*′=d*/dv，则当|v|充分小时，在v=0的某一侧存在Hopf分岔，且出现Hopf分岔时的极限环周期为

定理3是判断四维系统存在Hopf分岔的直接代数判据^[9]，应用该定理可避免过多烦琐的计算.由条件1所决定的曲面上各点就是发生Hopf分岔现象的分岔临界点.当系统在某工况下的参数满足条件1和2时，将出现周期振荡的失稳现象，即Hopf分岔现象，同时，该曲面也确定了系统的稳定域.从非线性动力系统结构稳定性分析的角度而言，当系统参数满足一定条件时，其结构稳定性会发生突变，出现系统结构失稳现象，从而导致系统表现出复杂非线性振荡行为.

对于四维系统，当发生双Hopf分岔时，其Jacobi矩阵可变换成如下形式：

其相应的特征方程：

应具有如下形式：

ω₁>0、ω₂>0为两对纯虚特征值的虚部，令式(4)和(5)中特征值λ对应的阶次系数相等，可得到四维非线性动力系统中发生双Hopf分岔时应满足的必要条件.

定理4：非线性系统y=f(x, v)中，x∈R⁴，v∈R¹为分岔参数，x=0为系统的平衡点，在平衡点时的系统Jacobi矩阵特征多项式为式(1)，若

成立，则当分岔参数v=0时，系统出现双Hopf分岔现象，且两个振荡频率ω₁和ω₂由下式获得：

当ω₁/ω₂为无理数时，双Hopf分岔称为非共振的，否则称为共振的，式(7)中，若p₂²-4p₀=0，则ω₁:ω₂=1:1，此时系统将出现1:1共振双Hopf分岔.

1.2 六维水轮机调节系统

对于四维以上的非线性系统，应采用下列定理.

定理5：非线性系统y=f(x, v)中，x∈Rⁿ，v∈R¹为分岔参数，x=0为系统的平衡点，在平衡点时的系统Jacobi矩阵特征多项式为

若当v=v₀时，有

1) Δ_i>0(i=1，2，…，n-2)，Δ_n-1=0，p₀(v₀)>0；

2) ω₀²UZ≠VW.

Δ_i(i=1，2，…，n)是由式(8)的系数所构造的Hurwitz行列式：

m=1, 2，…，n；若i>n，则p_n-1=0

n为偶数时：

n为奇数时:

当|v-v₀|充分小时，系统在v₀的某一侧存在Hopf分岔，定理6则进一步给出了判断高维非线性动力系统发生Hopf分岔的必要条件.

定理6：特征方程(8)有一对纯虚特征值，而其余n-2个特征值均具有负实部的充分必要条件是：

1) p_i>0(i=0，1，2，…，n-1)；

2) Δ_n-1=0，Δ_i>0(i=n-3，n-3，…).

定理5和6是判断非线性动力系统Hopf分岔存在性直接代数判据，通常用于维数较高的非线性系统.定理6条件2只要求一半Hurwitz行列式大于0，故其所涉及到的计算量小于定理5，因而实际应用更为简便.当满足定理6的两个条件时，非线性动力系统成为实际应用中常见的半单系统.由于分岔点是系统状态稳定的临界点，因此，分岔点的位置将决定系统在参数空间中的稳定域及稳定裕度.

对于六维系统，发生双Hopf分岔时，其Jacobi矩阵可变换为

其中，Q∈R².系统相应的特征方程:

应具有如下形式：

式中，ω₁>0，ω₂>0为两对纯虚特征值的虚部，方程λ²+aλ+b=0的两个根为矩阵Q的特征值.令式(10)、(11)中特征值λ的同阶次系数相等，可得到代数方程组：

若控制系统采用PID控制，则在系统特征方程系数p_i中，除ω₁、ω₂、a、b四个未知量外，还包括3个PID参数；应用中，根据先验知识或仿真结果，先确定其中一个控制参数，再通过方程组(12)计算出其余6个未知量.考虑到方程组(12)可能存在多解，实际计算中需将数值解代入到原方程组以获得稳态解，并计算系统Jacobi矩阵特征值进行验证.因此，式(12)可作为六维非线性系统双Hopf分岔的必要条件，用于水轮机调节系统六维非线性模型双Hopf分岔现象的研究.

2 数值计算与仿真 2.1 四维非线性系统模型

考虑发电机转子二阶动态模型和当前在水轮机调节系统中常用的PID控制，并结合文献[14, 15]中提出的刚性水击时的水轮机调节系统非线性模型，给出其四维非线性方程组：

其中：x、y、h、z为系统转速、接力器位置、水压、转速变化变量；T_w、T_a、T_y分别为水流惯性时间常数、机组惯性常数、接力器常数；e_xm、e_ym、e_hm、e_qxm、e_qhm、e_qym为线性水轮机组常数；e_g为电网自调节系数；r、m_g0分别为转速和力矩扰动; k_p、k_i、k_d为PID调节参数.

2.2 平衡点分析

上述方程组中，系统的平衡点位置与初始扰动状态有关，分3种情况：

1) 当r=0、m_g0=0时，系统有唯一的平衡点：(0，0，0，0)，系统无扰动；

2) 当r=0、m_g0≠0时，系统有唯一的平衡点：(0，m_g0/e_ym，0，m_g0/k_ie_ym)，存在负荷扰动；

3) 当r≠0、m_g0=0时，系统有唯一的平衡点：(r，(e_g-e_xm)r/e_ym，0，(e_g-e_xm)r/k_ie_ym)，存在频率扰动.

2.3 Hopf分岔

以某实际电站为例，某典型工况下，水轮机调节系统的参数分别为：T_w=2.0 s，T_a=6.87 s，e_xm=-1.067 3，e_g=0.4，e_ym=0.771 3，e_hm=1.717 9，e_qxm=-0.290 1，e_qhm=0.725 7，e_qym=0.818 4，T_y=0.1 s，则以上各系数为

用定理3来判断系统Hopf分岔存在性时，要求特征多项式系数p_i>0，由p₃>0，可得k_d < 6.364 6，即Hopf分岔区域在k_d < 6.364 6内.将式(13)代入式(2)，PID参数在出现Hopf分岔时应满足如下关系：

当以v_d=k_d-k_d^*或v_p=k_p-k_p^*作为分岔参数时，定理3中条件2左边为0，而右边不等于0，条件2成立.根据定理，系统将在|v|充分小时出现Hopf分岔现象，根据定理3和式(3)，发生Hopf分岔时极限环周期与PID控制参数的关系为

根据线性系统理论的Hurwitz稳定判据，可获得系统在PID参数空间中的稳定域，由式(14)给出系统在PID参数空间的分岔曲面，如图 2所示，该曲面将空间划分为稳定域和极限环两个性质不同的区域，曲面上每一点都可能是Hopf分岔临界点.Hopf分岔理论研究系统动力学特性的任务之一就是要确定这些分岔临界点，以确定系统稳定域边界及稳定域大小.为此，可在无扰动时由k_d取不同值在PI参数平面作出稳定域和极限环，其中，k_d分别取0、1、2和3，如图 3所示.

图 3中点1：k_d=2，k_i=1，由式(14)得k_p=4.424 7，v_p=0.图 4为机组转速x、接力器y和水头h在1点时的等幅低频振荡曲线, 由式(15)可计算出等幅振荡的周期为9.565 s, 与图 4中的仿真结果相一致.

图 5、6为图 4对应的系统变量相空间轨道，可以观察到在相空间一定区域中所形成的闭轨道极限环，系统呈扩散的增幅非线性振荡曲线，表明系统出现临界Hopf分岔.

2.4 Hopf分岔方向计算

上述图只给出了PID参数空间中Hopf分岔点的位置，但分岔方向及其稳定性的确定需要进一步计算.以v_p=k_p-k_p^*作分岔参数为例，在分岔点处系统Jacobi矩阵一对纯虚特征值为

将式(13)代入式(2)，可计算出系统纯虚特征值对v_p导数的实部为

同理，可计算出系统纯虚特征值对v_p导数虚部的表达式.在分岔点k_p=4.424 7、k_i=1.0、k_d=2.0处，可计算出a′(v₀)=0.138 5，再采用中心流形降阶方法^[9]可直接计算β=-2.075 2 < 0.根据定理2可判断系统Hopf分岔是超临界的，极限环将出现在分岔参数v₀>0一侧；而在v₀ < 0一侧系统是渐近稳定的.下面通过PI参数平面图 3上给出2、3点的参数对系统稳定性进行验证.

图 3中2点:k_d=2，k_i=1，取v_p=-0.4，则可得k_p=4.024 7，图 7为机组转速x、接力器y和水头h在2点时受2.5 Hz初始频率扰动下的振荡曲线，由式(15)可计算出振荡的周期为10.345 s，可知，系统在该组参数下逐渐趋于稳定，且振荡周期与计算值一致，如图 7所示.

图 3中3点：k_d =2，k_i=1，取v_p=-1，得到k_p=3.424 7, 图 8为机组转速x、接力器y和水头h在3点时受2.5 Hz初始频率扰动时的振荡曲线，由式(15)可计算出振荡的周期为10.986 s，系统在该组参数下经过5个振荡周期后能迅速收敛于稳定点，该参数下系统是稳定的.

图 4~8表明了系统的分岔方向及稳定性.当水轮机调节系统中采用PID(或增益调度PID)策略时，应将PID参数设置在图 3中极限环曲线下方且远离曲线的区域，否则系统中将有可能出现持续的或不稳定的增幅非线性振荡，而无法获得稳定的控制效果.

2.5 双Hopf分岔

当k_d=6.364 6，且k_p与k_i关系满足式(6)时，系统会出现双Hopf分岔现象.式(6)中，若取k_i=3.0，则k_p=4.680 3，PID参数空间中双Hopf分岔临界点之一为k_p=4.680 0、k_i=3.0、k_d=6.364 6.计算得到此时系统Jacobi矩阵的两对特征值分别为-2.748 0^-13±i0.706 1和-1.332 7^-14±i2.157 3，考虑到数值计算中的允许误差，可认为两对共轭特征值的实部为0，系统存在两对纯虚特征根.图 9给出两对复共轭特征值随k_p的变化在双Hopf分岔临界点处同时穿越虚轴的情况，两对复共轭特征值穿越方向相异，表明发生双Hopf分岔前后系统处于不稳定状态，与图 1(b)中的情况对应.

图 10为系统出现双Hopf分岔现象状态变量的时间振荡曲线，从中可观察到，从平衡点分岔出来的将不再是周期轨道，而可能是拟(准)周期轨道，两个不同频率的振荡相互叠加形成更为复杂的非线性振荡波形，该波形与非线性系统的Neimark-Sacker分岔(环面分岔)中的调制波形非常相似.

当发生双Hopf分岔时，在系统状态空间中将会出现闭轨道双极限环以及轮胎面或指环形的吸引域，如图 11所示，包含丰富的动力学特征，是一种比相空间中极限环更为复杂的非线性动态行为.

2.6 六维非线性系统模型双Hopf分岔

当取弹性水击模型时，上述四维可扩展到六维^{[14, 15]}，水力系统参数T_r=1.35，h_w=2，方程组为

水轮机参数取值同上，且平衡点分析、Hopf分岔、分岔方向及稳定性可依据文中给出的公式和定理进行计算，分析过程及结果与四维系统的类似.根据六维非线性系统出现双Hopf分岔的必要条件式(12)，给定k_i时计算发生双Hopf分岔的a、b、k_p、k_d、ω₁和ω₂等值.表 1~3分别给出了无扰动、2.5 Hz频率扰动和30%负荷扰动时双Hopf分岔临界点及系统特征值的计算结果，在一定的PID参数范围内，系统Jacobi矩阵出现两对纯虚共轭特征值，且两对纯虚部值是系统发生双Hopf分岔时的两个振荡频率ω₁和ω₂，当系统以两个不同频率进行振荡时，将呈现出更为复杂的动态特性.

图 12、13为六维系统在PID控制参数分别为k_p=1.315、k_i=1、k_d=5.591并受到2.5 Hz频率扰动的初始条件下，双Hopf分岔两个不同频率的状态变量振荡曲线相互迭加的时域响应图和对应的相空间吸引域.在流体力学关于湍流发生机制的理论中，这种拟(准)周期轨道所引发的轮胎面吸引域事实上也是发生湍流Landau-Hopf^[16]理论分析的出发点.

3 结论

本文基于非线性动力系统的分岔理论，分别对水轮机调节系统的刚性和弹性水击连续时间非线性解析模型进行了理论分析，给出了系统在PID参数空间中不同情况下的稳定域范围，证明了所研究的非线性系统均存在超临界Hopf分岔现象，由于其经分岔产生的极限环是稳定的，因而可进一步观察系统的分岔特性，包括分岔参数的临界值、在稳定域内的平衡点位置以及发生分岔时极限环的方向等；给出了判断系统出现双Hopf分岔的必要条件，对水轮机调节系统的四维和六维非线性模型进行了双Hopf分岔的存在性的计算，对于考虑水轮机的传递系数随工况变化的非线性动态模型，无论是刚性水击还是弹性水击模型，都存在双Hopf分岔现象，通过计算给出了发生双Hopf分岔时的参数值.最后, 通过数值计算对分岔点两侧系统稳定性的变化情况进行仿真，验证了理论分析结果.文中的分析方法适用于一般水力发电系统的稳定运行分析和研究, 从而为水轮机调节系统中控制参数的正确选择提供了一条新的理论分析途径.