CN114329835A - 一种多机共管道水轮机调节系统稳定域确定方法及系统 - Google Patents

Abstract

本发明提供一种多机共管道水轮机调节系统稳定域确定方法及系统，包括：确定高阶状态空间模型；模型通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；确定系统矩阵的特征值；当特征值的实部均小于零时，水轮机调节系统处于稳定状态；通过变步长和二分法确定水轮机调节系统的稳定域；采用变步长方式确定稳定域的第一种稳定边界，第一种稳定边界为三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定稳定域的第二种稳定边界，第二种稳定边界为三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。本发明实现对系统稳定域的快速计算。

Description

一种多机共管道水轮机调节系统稳定域确定方法及系统

技术领域

本发明属于水轮机调节系统建模与稳定性分析技术领域，更具体地，涉及一种多机共管道水轮机调节系统稳定域确定方法及系统。

背景技术

过去几年，中国以水电为主可再生能源发电量持续增加，并规划建设雅鲁藏布江下游水电能源基地和金沙江上下大型清洁能源基地等。大力发展可再生能源与提高间歇性能源消纳能力已经成为全球能源改革趋势，水电能源在其中起着关键作用。水电能源作为运用规模最大的可再生能源，在许多电力系统中承担着调节、平衡与电能消纳的任务。在许多水能丰富的地区如中国的西南，更是形成了以水电为主导的电网结构。扩大水电生产规模与提高水电运行质量是支持长远发展可再生能源战略的结构基础。如何确保水电系统灵活持续稳定性是水力发电中最为基础与重要的工作。

水轮机调节系统是一个具有非最小相位特征水-机-电耦合的复杂系统，多机共管道水电系统相比单机系统水力耦合更为复杂，运行方式也更为灵活，其不仅可能运行在不同的水轮机工况，还可能以不同的机组组合和不同的机组运行。原本整定的PI控制参数在复杂工况下难以满足大工况变化的稳定性要求，极端情况甚至导致系统失稳。因此对于水电系统需要明确全水轮机工况下稳定性，对于多机公管道水轮机调节系统还应考虑灵活运行情景条件下的稳定性。

目前，水轮机调节系统稳定性分析方法主要源自经典控制理论或现代控制理论。常用的水轮机调节系统稳定性分析方法有Routh-Hurwitz准则，根轨迹，Hopf分叉理论以及Lyapunov理论等。但计算系统稳定域主要采用Hurwitz准则，这种方法计需要进行一系列的规划，对于低阶模型较为适用。但是如果模型阶数较高，求解难度将显著增加，计算速度将显著变慢，对于过高阶数的模型还可能无法求解。这严重制约了水轮机调节系统灵活运行条件下的稳定性研究。深入研究多机系统多水轮机工况和多种运行情景的条件的稳定性，充分发挥水电系统灵活运行的优势，亟需解决高阶系统稳定域计算问题与批量求解系统稳定域耗时较长的问题。因此需要一种迫切需要一种可以灵活，快速计算多机水力耦合水轮机调节系统稳定域的方法。

发明内容

针对现有技术的缺陷，本发明的目的在于提供一种多机共管道水轮机调节系统稳定域确定方法及系统，旨在解决现有解决高阶系统稳定域计算问题与批量求解系统稳定域耗时较长的问题。

为实现上述目的，本发明提供了一种多机共管道水轮机调节系统稳定域确定方法，包括如下步骤：

确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机组四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

在一个可选的示例中，所述PID型调速器包括PID控制器；

所述PID控制器的微分方程表示为：

式中，K_p为比例系数，K_i为积分系数，K_d为微分系数，T_D为微分时间常数，b_p为永态转差系数，x_p,x_i,x_d分别为比例信号、积分信号、微分信号，x为机组转速偏差相对值，c为设定转速相对值，y₁为辅助接力器输出信号；

辅助接力器和主接力器的微分方程为：

式中，T_y1辅助接力器反应时间常数，T_y为主接力器反应时间常数，y为导叶开度偏差相对值，s为拉普拉斯算子。

在一个可选的示例中，所述引水系统包括：压力管道、尾水管、调压室、汇合管道以及尾水隧洞；

其中压力管道、尾水管以及尾水隧洞的微分方程分别为：

压力管道微分方程：

尾水管微分方程：

尾水隧洞微分方程：

式中，h_b1,h_x1,h_p1,q₂₁均为引入计算二阶近似弹性水击而引入的计算中间状态变量，T_r1,T_r2,T_r3分别为压力管道、尾水管和尾水隧洞的水击相长；h_w1，h_w2，h_w3分别为压力管道、尾水管和尾水隧洞的管道特征系数；h_b为压力管道末端断面的水头偏差相对值，h_x为尾水管入口水头的偏差相对值，q_t为水轮机流量的偏差相对值，q₂为尾水管出口流量的偏差相对值，h_p为尾水隧洞入口流量的偏差相对值；a₃为尾水隧洞的管道损失常数；

汇合管道与调压室微分方程为：

调压室微分方程：

汇合管道微分方程：

式中，q_c为调压室流量的偏差相对值，h_p为调压室水头的偏差相对值，T_c为调压室水头时间常数，T_q为调压室流量时间常数；h_i表示第i根管道在汇合断面的水头，i∈(1,2,…,n,n+1)，q_i表示第i根管道在汇合断面的流量，n+1管道是共用的总管道。

在一个可选的示例中，所述水轮机的模型表示为：

式中，

m_t为水轮机力矩偏差相对值；q_t为水轮机流量偏差相对值；x为转速偏差相对值，h_t为水轮机水头偏差相对值；y为导叶开度偏差相对值；e_y为水轮机力矩对导叶开度传递系数；e_x为水轮机力矩对转速传递系数；e_h为水轮机力矩对工作水头传递系数；e_qy为水轮机流量对导叶开度传递系数；e_qx为水轮机流量对转速传递系数；e_qh为水轮机流量对工作水头传递系数。

在一个可选的示例中，所述发电机组的微分方程表示为：

式中，T_ae为等效机组惯性时间常数，m_te为机组等效主动力矩的偏差相对值，m_g0为等效负载力矩偏差相对值，e_ne为等效发电机负载自调节系数；k_i为第i台机组的额定容量占比，T_ai为第i台机组的惯性时间常数，e_g为发电机负载自调节系数，e_xi为第i台机组力矩对转速传递系数，m_ti为第i台机组的主动力矩偏差相对值，n表示一共有n台机组。

在一个可选的示例中，当所述机组为两台时，所述高阶状态空间模型为：

其中，x＝[x,y₁,y'₁,y,y',x_i,x'_i,x_d,x'_d,h_x1,h_x,h'_x1,h'_x,q₂₁,q₂,q'₂₁,q'₂,h_b1,h_b,h'_b1,h'_b,h_p1,h_p,q_c]为状态变量，式中，没有上标的变量为第一台机组单元的变量，有上标“'”的变量表示第二台机组单元的变量；

a_1,1＝-e_ne/T_ae,a_1,4＝ke_y/T_ae,a_1,5＝k'e_y'/T_ae,a_1,11＝-ke_h/T_ae；

a_1,13＝-k'e'_h/T_ae,a_1,19＝ke_h/T_ae,a_1,21＝k'e'_h；

a_2,1＝-K_p/T_y1,a_2,2＝-(1+K_pb_p)/T_y1,a_2,6＝1/T_y1,a_2,8＝1/T_y1；

a_3,1＝-K'_p/T'_y1,a_3,3＝-(1+K'_pb'_p)/T'_y1,a_3,9＝1/T'_y1,a_3,12＝1/T'_y1a_4,2＝1/T_y,a_4,4＝-1/T_y；

a_5,3＝1/T'_y,a_5,5＝-1/T'_y；

a_6,1＝-K_i,a_6,2＝-K_ib_p；a_7,1＝-K'_i,a_7,3＝-K'_ib'_p；

a_8,1＝(K_de_ne)/(T_dT_ae),a_8,4＝-(kK_de_y)/(T_dT_ae),a_8,5＝-(k'K_de'_y)/(T_dT_ae),a_8,8＝-1/T_d；

a_8,11＝(kK_de_h)/(T_dT_ae),a_8,13＝(k'K_de'_h)/(T_dT_ae),a_8,19＝(k K_de_h)/(T_dT_ae),a_8,21＝-(k'K_de'_h)/(T_dT_ae)；

a_9,1＝(K'_de_ne)/(T'_dT_ae),a_9,4＝-(kK'_de_y)/(T'_dT_ae),a_9,5＝-(k'K'_de'_y)/(T'_dT_ae),a_9,9＝-1/T'_d；

a_9,11＝(kK'_de_h)/(T'_dT_ae),a_9,13＝(k'K'_de'_h)/(T'_dT_ae),a_9,19＝-(kK'_de_h)/(T'_dT_ae),a_9,21＝(k'K'_de'_h)/(T'_dT_ae)；

a_11,1＝8h_w2e_qx/T_r2,a_11,4＝8h_w2e_qy/T_r2,a_11,10＝1,a_11,11＝-8h_w2e_qh/T_r2,a_11,19＝8h_w2e_qh/T_r2；

a_13,1＝8h'_w2e'_qx/T_r'₂,a_13,5＝8h'_w2e'_qy/T_r2,a_13,12＝1,a_13,13＝-8h'_w2e'_qh/T'_r2,a_13,21＝8h'_w2e'_qh/T'_r2；

a_15,14＝1,a_15,23＝-2/(h_w2T_r2)；

a_17，16＝1，a_17，23＝-2/(h′_w2T′_r2)；

a_19,1＝-8h_w1e_qx/T_r1,a_19,4＝-8h_w1e_qy/T_r1,a_19,11＝8h_w1e_qh/T_r1,a_19,18＝1,a_19,19＝-8h_w1e_qh/T_r1；

a_21,1＝-8h'_w1e'_qx/T'_r1,a_21,5＝-8h'_w1e'_qy/T'_r1,a_21,13＝8h'_w1e'_qh/T'_r1,a_21,20＝1,a_21,21＝-8h'_w1e'_qh/T'_r1；

a_23,21＝8h_w3/T_r3,a_23,17＝8h_w3/T_r3,a_23,22＝1,a_23,23＝-4a₃/T_r3,a_23,24＝-8h_w3/T_r3；

a_24,15＝8h_w3T_c/(T_r3T_q),a_24,17＝8h_w3T_c/(T_r3T_q),a_24,22＝T_c/T_q；

a_24,23＝-4a₃T_c/(T_r3T_q),a_24,24＝-(8h_w3T_c+T_r3)/(T_r3T_q)；

B为输入矩阵，所述系统矩阵为：

在一个可选的示例中，通过变步长和二分法确定所述水轮机调节系统的稳定域，具体为：

当微分系数K_d为预设固定值时，将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中：

Step 1：采用变步长的方法循环选取比例系数K_p和积分系数K_i的取值为(0，K_{i_0})带入系统矩阵中计算其特征值以判断系统稳定性，进而求得PI临界稳定线的起点(0，K_{i_0})；其中，坐标原点(0，0)到(0，K_{i_0})之间的PI临界稳定线为第一种稳定边界，K_{i_0}为积分系数K_i在第一种稳定边界的最大值；

Step 2：按照固定的ΔK_p步长，利用二分法在所述临界稳定线上的上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，具体步骤为：

Step2.1：给定计算精度e,K_p搜索步长dK_p＝0.01，搜索邻域的角度θ范围为初始范围[l,u]初始l＝-pi/2，u＝pi/2；K_i的基准值是K_{i_0}；

Step2.2：在[l，u]区间取中点θ₁；

Step2.3：在θ与dK_p构成三角函数关系：K_i＝dK_p*tanθ₁，计算并判断(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)对应的系统状态；

Step2.4：若系统稳定将θ₁作为θ区间的下限即l＝θ₁，若系统不稳定则将θ₁作为θ区间的上限即u＝θ₁，重复Step2.2～Step2.4，直到|u-l|<e；

Step2.5：将Step2.4满足精度条件的PI参数即(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)作为新Step2.1的基准点，重复Step2.1～Step2.5连续求取新的临界稳定参数，直到K_i<0停止；

Step 3：采用变步长的方法循环选取的(K_p，0)带入系统矩阵中计算特征值以判断系统稳定性进而，求得PI临界稳定线的终点(K_{p_end}，0)；其中，坐标原点(0，0)到(K_{p_end}，0)之间的PI临界稳定线为第一种稳定边界，K_{p_end}为比例系数K_p在第一种稳定边界的最大值；

将Step 1～Step 3求取的PI临界稳定线的数据结合起来，其内部的区域属于所求的稳定域；

将比例系数K_p和积分系数K_i中任一系数设为预设固定值，参见上述Step 1～Step3的方法求取微分系数K_d的上限；

改变微分系数K_d的预设固定值，参见上述Step 1～Step 3求取不同微分系数K_d下的稳定域，将其组合得到系统的稳定域。

在一个可选的示例中，所述步骤Step 1具体包括如下步骤：

Step 1.1：以坐标原点(0，0)为基准点，给定初始步长dk_i＝0.1，在纵轴上找到下一个点判断系统状态，如果稳定则将系统的将步长dk_i＝5*dk_i，并以该点纵坐标为基准更新下一个纵坐标；如果系统不稳定，则将系统步长变为dk_i＝0.2*dk_i，以上一个稳定的点为基础更新纵坐标；

Step 1.2：重复Step 1.1的更新步骤，直到两次纵坐标的距离小于设定的精度，此时的纵坐标值即为PI临界稳定线的起点(0，K_{i_0})；

所述步骤Step 3具体包括如下步骤：

Step 3.1：以坐标原点(0，0)为基准点，给定初始步长dk_p＝0.1，在横轴上找到下一个点判断系统状态，如果稳定则将系统的将步长dk_p＝5*dk_p，并以该点横坐标为基准更新下一个横坐标；如果系统不稳定，则将系统步长变为dk_p＝0.2*dk_p，以上一个稳定的点为基础更新横坐标；

Step 3.2：重复Step 3.1的更新步骤，直到两次横坐标的距离小于设定的精度，此时的横坐标值即为PI临界稳定线的终点(K_{p_end}，0)。

第二方面，本发明提供了一种多机共管道水轮机调节系统稳定域确定系统，包括：

系统模型确定单元，用于确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

特征值确定单元，用于确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

稳定域确定单元，用于通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

在一个可选的示例中，所述稳定域确定单元通过变步长和二分法确定所述水轮机调节系统的稳定域，具体为：

当微分系数K_d为预设固定值时，将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中：

Step 1：采用变步长的方法循环选取比例系数K_p和积分系数K_i的取值为(0，K_{i_0})带入系统矩阵中计算其特征值以判断系统稳定性，进而求得PI临界稳定线的起点(0，K_{i_0})；其中，坐标原点(0，0)到(0，K_{i_0})之间的PI临界稳定线为第一种稳定边界，K_{i_0}为积分系数K_i在第一种稳定边界的最大值；

Step 2：按照固定的ΔK_p步长，利用二分法在所述临界稳定线上的上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，具体步骤为：

Step2.1：给定计算精度e,K_p搜索步长dK_p＝0.01，搜索邻域的角度θ范围为初始范围[l,u]初始l＝-pi/2，u＝pi/2；K_i的基准值是K_{i_0}；

Step2.2：在[l，u]区间取中点θ₁；

Step2.3：在θ与dK_p构成三角函数关系：K_i＝dK_p*tanθ₁，计算并判断(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)对应的系统状态；

Step2.4：若系统稳定将θ₁作为θ区间的下限即l＝θ₁，若系统不稳定则将θ₁作为θ区间的上限即u＝θ₁，重复Step2.2～Step2.4，直到|u-l|<e；

Step2.5：将Step2.4满足精度条件的PI参数即(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)作为新Step2.1的基准点，重复Step2.1～Step2.5连续求取新的临界稳定参数，直到K_i<0停止；

Step 3：采用变步长的方法循环选取的(K_p，0)带入系统矩阵中计算特征值以判断系统稳定性进而，求得PI临界稳定线的终点(K_{p_end}，0)；其中，坐标原点(0，0)到(K_{p_end}，0)之间的PI临界稳定线为第一种稳定边界，K_{p_end}为比例系数K_p在第一种稳定边界的最大值；

将Step 1～Step 3求取的PI临界稳定线的数据结合起来，其内部的区域属于所求的稳定域；

将比例系数K_p和积分系数K_i中任一系数设为预设固定值，参见上述Step 1～Step3的方法求取微分系数K_d的上限；

改变微分系数K_d的预设固定值，参见上述Step 1～Step 3求取不同微分系数K_d下的稳定域，将其组合得到系统的稳定域。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，具有以下有益效果：

1)本发明提供一种多机共管道水轮机调节系统稳定域确定方法及系统，建立了双机共尾水隧洞水轮机调节系统的高阶状态空间模型，便于分析水电机组在在频率控制模式下小扰动稳定性和动态响应特性。

2)本发明提供一种多机共管道水轮机调节系统稳定域确定方法及系统，根据状态矩阵计算系统特征值，并依全部特征值实部是否均为负数可以直接判断特征值与系统稳定性的关系。

3)本发明提供一种多机共管道水轮机调节系统稳定域确定方法及系统，将调速器的PI控制参数表示为直角坐标系中的点，并依据变步长搜索临界稳定线起点(0，K_{i_0})和终点(K_{p_end}，0)以及二分法连续邻域搜索临界稳定点(K_p，K_i)，可实现精度高且计算速度快的高阶系统稳定域求解。

附图说明

图1是本发明实施例提供的多机共管道水轮机调节系统稳定域确定方法流程图；

图2是本发明实施例提供的高阶水轮机调节系统模型框图；

图3是本发明实施例提供的多机共管道水轮机调节系统稳定域求解的步骤示意图；

图4是本发明实施例提供的系统稳定域搜索策略图；

图5是本发明实施例计算的系统稳定域结果示意图；

图6是本发明实施例提供的多机共管道水轮机调节系统稳定域确定系统架构图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

本发明涉及一种多机共管道水轮机调节系统稳定域快速计算方法，属于水轮机调节系统建模与稳定性分析技术领域。该方法具体操作步骤为：建立多机共管道水轮机调节系统状态空间模型，并依据系统的状态矩阵所有的特征值的实部是否均为负数判断对应系统是否稳定；将调速器临界稳定的PI控制参数作为求解目标，并将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中。首先采用变步长的方法循环选取的(0，K_i)带入状态矩阵中计算特征值判断系统稳定性进而求得PI临界稳定线的起点(0，K_{i_0})；然后按照固定的ΔK_p步长，利用二分法在上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，并以这种形式连续求解PI稳定域的稳定边界；最后利用变步长的方法循环选取的(K_p，0)带入状态矩阵中计算特征值判断系统稳定性进而求取求得PI临界稳定线的终点(K_{p_end},0)；结合起点、终点和连续搜索结果得到调速器PI参数的稳定域。本发明在保证高计算精度的同时显著的提升了计算速度，可以为研究多机共管道水轮机调节系统在灵活运行的条件下的稳定性提供有力的计算方法。

本发明的目的在于提供一种快速，灵活的多机水轮机调节系统控制参数稳定域计算方法，具有计算速度快、精度较高的优点。

图1是本发明实施例提供的多机共管道水轮机调节系统稳定域确定方法流程图；如图1所示，包括如下步骤：

S101，确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机组四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

S102，确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

S103，通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

具体如下实施例所述：本实例以我国某水电站的双机共尾水水电系统为研究对象，具体系统参数给出如表1所示。仿真工况为机组以额定出力并网运行时发生+0.1p.u.的负荷扰动。

表1水电系统基本参数

本发明涉及一种多机共管道水轮机调节系统稳定域快速计算方法，包括以下步骤：

(1)建立包含PID型调速器、引水系统、水轮机、发电机的双机共尾水水轮机调节系统高阶状态空间模型，结构如附图2所示。

1)PID型调速器模型

PID型调速器包含PID控制器和和液压放大部分。其中PID控制器的微分方程表示为：

式中，K_p为比例系数，K_i为积分系数，K_d为微分系数，T_D为微分时间常数，b_p为永态转差系数，x_p,x_i,x_d分别为比例信号、积分信号、微分信号，x为机组转速偏差相对值，c为设定转速相对值，y₁为辅助接力器输出信号；

辅助接力器和主接力器的微分方程为：

式中，T_y1辅助接力器反应时间常数，T_y为主接力器反应时间常数，y为导叶开度偏差相对值。

2)引水系统模型

多机共管道的水系统包含压力管道，尾水管，调压室，汇合管道，以及尾水隧洞。其中压力管道，尾水管以及尾水隧洞分别为：

压力管道

尾水管

尾水隧洞

式中，h_b1,h_x1,h_p1,q₂₁引入计算二阶近似弹性水击而引入的计算中间状态变量，T_r1,T_r2,T_r3分别为压力管道、尾水管和尾水隧洞的水击相长；h_w1，h_w2，h_w3分别为压力管道、尾水管和尾水隧洞的管道特征系数；h_b为压力管道末端断面的水头偏差相对值，h_x为尾水管入口水头的偏差相对值，q_t为水轮机流量的偏差相对值，q₂为尾水管出口流量的偏差相对值，h_p为尾水隧洞入口流量的偏差相对值；a₃为尾水隧洞的管道损失常数。

汇合管道与调压室的微分方程为：

调压室

汇合管道

式中，q_c为调压室流量的偏差相对值，h_p为调压室水头的偏差相对值，Tc为调压室水头时间常数，T_q为调压室流量时间常数；h_i表示第i根管道在汇合断面的水头，q_i表示第i根管道在汇合断面的流量，n+1管道是共用的总管道。

3)水轮机模型

水轮机用IEEE六参数模型表示为：

式中，

m_t为水轮机力矩偏差相对值；e_y为水轮机力矩对导叶开度传递系数；e_x为水轮机力矩对转速传递系数；e_h为水轮机力矩对工作水头传递系数；e_qy为水轮机流量对导叶开度传递系数；e_qx为水轮机流量对转速传递系数；e_qh为水轮机流量对工作水头传递系数。

4)发电机负载模型

多台发电机及负荷之间近似认为刚性联结，采用等效的一阶转子运动方程，微分方程表示为：

式中，T_ae为等效机组惯性时间常数，m_te为机组等效主动力矩的偏差相对值，m_g0为等效负载力矩偏差相对值，e_ne为等效发电机负载自调节系数；k_i为第i台机组的额定容量占比，T_ai为第i台机组的惯性时间常数，e_g为发电机负载自调节系数，e_xi为第i台机组力矩对转速传递系数，m_ti为第i台机组的主动力矩偏差相对值，n表示一共有n台机组。

5)多机共管道水轮机调节系统高阶状态空间模型

以机组有2台为例进行举例说明：选取

x＝[x,y₁,y'₁,y,y',x_i,x'_i,x_d,x'_d,h_x1,h_x,h'_x1,h'_x,q₂₁,q₂,q'₂₁,q'₂,h_b1,h_b,h'_b1,h'_b,h_p1,h_p,q_c]为状态变量，式中，没有上标的变量为第一台机组单元的变量，有上标“'”的变量表示第二台机组单元的变量；

则多机共管道水轮机调节系统高阶状态空间模型表示为：

a_1,1＝-e_ne/T_ae,a_1,4＝ke_y/T_ae,a_1,5＝k'e_y'/T_ae,a_1,11＝-ke_h/T_ae；

a_1,13＝-k'e'_h/T_ae,a_1,19＝ke_h/T_ae,a_1,21＝k'e'_h；

a_2,1＝-K_p/T_y1,a_2,2＝-(1+K_pb_p)/T_y1,a_2,6＝1/T_y1,a_2,8＝1/T_y1；

a_3,1＝-K'_p/T'_y1,a_3,3＝-(1+K'_pb'_p)/T'_y1,a_3,9＝1/T'_y1,a_3,12＝1/T'_y1a_4,2＝1/T_y,a_4,4＝-1/T_y；

a_5,3＝1/T'_y,a_5,5＝-1/T'_y；

a_6,1＝-K_i,a_6,2＝-K_ib_p；a_7,1＝-K'_i,a_7,3＝-K'_ib'_p；

a_8,1＝(K_de_ne)/(T_dT_ae),a_8,4＝-(kK_de_y)/(T_dT_ae),a_8,5＝-(k'K_de'_y)/(T_dT_ae),a_8,8＝-1/T_d；

a_8,11＝(kK_de_h)/(T_dT_ae),a_8,13＝(k'K_de'_h)/(T_dT_ae),a_8,19＝(k K_de_h)/(T_dT_ae),a_8,21＝-(k'K_de'_h)/(T_dT_ae)；

a_9,1＝(K'_de_ne)/(T'_dT_ae),a_9,4＝-(kK'_de_y)/(T'_dT_ae),a_9,5＝-(k'K'_de'_y)/(T_d'T_ae),a_9,9＝-1/T'_d；

a_9,11＝(kK'_de_h)/(T'_dT_ae),a_9,13＝(k'K'_de'_h)/(T'_dT_ae),a_9,19＝-(kK'_de_h)/(T'_dT_ae),a_9,21＝(k'K'_de'_h)/(T'_dT_ae)；

a_11,1＝8h_w2e_qx/T_r2,a_11,4＝8h_w2e_qy/T_r2,a_11,10＝1,a_11,11＝-8h_w2e_qh/T_r2,a_11,19＝8h_w2e_qh/T_r2；

a_13,1＝8h'_w2e'_qx/T'_r2,a_13,5＝8h'_w2e'_qy/T_r2,a_13,12＝1,a_13,13＝-8h'_w2e'_qh/T'_r2,a_13,21＝8h'_w2e'_qh/T'_r2；

a_15,14＝1,a_15,23＝-2/(h_w2T_r2)；

a_17，16＝1，a_17，23＝-2/(h′_w2T′_r2)；

a_19,1＝-8h_w1e_qx/T_r1,a_19,4＝-8h_w1e_qy/T_r1,a_19,11＝8h_w1e_qh/T_r1,a_19,18＝1,a_19,19＝-8h_w1e_qh/T_r1；

a_21,1＝-8h'_w1e'_qx/T'_r1,a_21,5＝-8h'_w1e'_qy/T'_r1,a_21,13＝8h'_w1e'_qh/T'_r1,a_21,20＝1,a_21,21＝-8h'_w1e'_qh/T'_r1；

a_23,21＝8h_w3/T_r3,a_23,17＝8h_w3/T_r3,a_23,22＝1,a_23,23＝-4a₃/T_r3,a_23,24＝-8h_w3/T_r3；

a_24,15＝8h_w3T_c/(T_r3T_q),a_24,17＝8h_w3T_c/(T_r3T_q),a_24,22＝T_c/T_q；

a_24,23＝-4a₃T_c/(T_r3T_q),a_24,24＝-(8h_w3T_c+T_r3)/(T_r3T_q)；

B为输入矩阵，对本稳定域求解方法没有影响，故不将系数一一列出。

此时，系统矩阵A为：

(2)根据系统状态矩阵的所有特征值实部的是否均为负号判断系统是否稳定。

利用状态方程的系统矩阵A可计算系统矩阵的特征值为：

λ_i＝α_i+jβ_i i＝1,2,…,24

式中，λ_i为系统矩阵的第i个特征值，α_i表示系统矩阵第i个特征值的实部，β_i表示系统矩阵第i个特征值的虚部。

系统稳定的充分条件是所有α_i<0 i＝1,2,…,24；如果存在任意的一个α_i>0 i＝1,2,…,24则系统不稳定。

(3)将求取临界稳定的PI控制参数的作为求取目标，并将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中。

(4)参见图3，首先，Step 1采用变步长的方法循环选取的(0，K_i)带入状态矩阵中计算特征值判断系统稳定性进而求得PI临界稳定线的起点(0，K_{i_0})。其次，Step 2按照固定的ΔK_p步长，利用二分法在上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，并以这种形式连续求解PI稳定域的稳定边界。Step 3利用变步长的方法循环选取的(K_p，0)带入状态矩阵中计算特征值判断系统稳定性进而求取求得PI临界稳定线的终点(K_{p_end},0)。最后，结合起点，终点和连续邻域搜索的结果得到调速器PI参数的稳定域。具体流程图如附图4所示。

1)采用变步长的方法循环选取的(0，K_i)带入状态矩阵中计算特征值判断系统稳定性进而求得PI临界稳定线的起点(0，K_{i_0})，其步骤为：

Step 1.1：坐标原点(0，0)为基准点，给定初始步长dk_i＝0.1，在纵轴上找到下一个点判断系统状态，如果稳定则将系统的将步长dk_i＝5*dk_i，并以该点纵坐标为基准更新下一个纵坐标；如果系统不稳定，则将系统步长变为dk_i＝0.2*dk_i，以上一个稳定的点为基础更新纵坐标。

Step 1.2：重复Step 1.1的更新步骤，直到两次纵坐标的距离小于设定的精度，此时的纵坐标值即为K_p-K_i临界稳定线的起点(0，K_{i_0})。

2)按照固定的ΔK_p步长，利用二分法在上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，并以这种形式连续求解PI稳定域的稳定边界，其步骤为：

Step2.1：给定计算精度e,K_p搜索步长dK_p＝0.01，搜索邻域的角度θ范围为初始范围[l,u]初始l＝-pi/2，u＝pi/2；K_i的基准值k_i＝K_{i_0}；初始的K_i基准值为Step 1输出的纵坐标。

Step2.2：在[l，u]区间取中点θ₁，

Step2.3：在θ与dK_p构成三角函数关系：K_i＝dK_p*tanθ₁，计算并判断(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)对应的系统状态。

Step2.4：若系统稳定将θ₁作为θ区间的下限即l＝θ₁，若系统不稳定则将θ₁作为θ区间的上限即u＝θ₁，重复Step2.2～Step2.4，直到|u-l|<e。

Step2.5：将Step2.4满足精度条件的PI参数即(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)作为新Step2.1的基准点，重复Step2.1～Step2.5连续求取新的临界稳定参数，直到K_i<0停止。

3)利用变步长的方法循环选取的(K_p，0)带入状态矩阵中计算特征值判断系统稳定性进而求取求得PI临界稳定线的终点(K_{p_end},0)，其步骤为Step 3.1：坐标原点(0，0)为基准点，给定初始步长dk_p＝0.1，在横轴上找到下一个点判断系统状态，如果稳定则将系统的步长变为dk_p＝5*dk_p，并以该点横坐标为基准更新下一个横坐标；如果系统不稳定，则将系统步长变为dk_p＝0.2*dk_p，以上一个稳定的点为基础更新横坐标。

Step 3.2：重复Step 3.1的更新步骤，直到两次横坐标的距离小于设定的精度，此时的横坐标值即为K_p-K_i临界稳定线的终点(K_{p_end}，0)。

4)将Step 1～Step 3的数据结合起来，其内部的区域为所求的稳定域，参见图5所示。

图6是本发明实施例提供的多机共管道水轮机调节系统稳定域确定系统架构图，如图6所示，包括：

系统模型确定单元610，用于确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机组四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

特征值确定单元620，用于确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

稳定域确定单元630，用于通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

可以理解的是，图6中各个单元的功能实现可参见前述方法实施例中的详细介绍，在此不做赘述。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种多机共管道水轮机调节系统稳定域确定方法，其特征在于，包括如下步骤：

确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机组四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

2.根据权利要求1所述的稳定域确定方法，其特征在于，所述PID型调速器包括PID控制器；

所述PID控制器的微分方程表示为：

式中，K_p为比例系数，K_i为积分系数，K_d为微分系数，T_D为微分时间常数，b_p为永态转差系数，x_p,x_i,x_d分别为比例信号、积分信号、微分信号，x为机组转速偏差相对值，c为设定转速相对值，y₁为辅助接力器输出信号；

辅助接力器和主接力器的微分方程为：

式中，T_y1辅助接力器反应时间常数，T_y为主接力器反应时间常数，y为导叶开度偏差相对值，s为拉普拉斯算子。

3.根据权利要求2所述的稳定域确定方法，其特征在于，所述引水系统包括：压力管道、尾水管、调压室、汇合管道以及尾水隧洞；

其中压力管道、尾水管以及尾水隧洞的微分方程分别为：

压力管道微分方程：

尾水管微分方程：

尾水隧洞微分方程：

式中，h_b1,h_x1,h_p1,q₂₁均为引入计算二阶近似弹性水击而引入的计算中间状态变量，T_r1,T_r2,T_r3分别为压力管道、尾水管和尾水隧洞的水击相长；h_w1，h_w2，h_w3分别为压力管道、尾水管和尾水隧洞的管道特征系数；h_b为压力管道末端断面的水头偏差相对值，h_x为尾水管入口水头的偏差相对值，q_t为水轮机流量的偏差相对值，q₂为尾水管出口流量的偏差相对值，h_p为尾水隧洞入口流量的偏差相对值；a₃为尾水隧洞的管道损失常数；

汇合管道与调压室微分方程为：

调压室微分方程：

汇合管道微分方程：

式中，q_c为调压室流量的偏差相对值，h_p为调压室水头的偏差相对值，T_c为调压室水头时间常数，T_q为调压室流量时间常数；h_i表示第i根管道在汇合断面的水头，i∈(1,2,…,n,n+1)，q_i表示第i根管道在汇合断面的流量，n+1管道是共用的总管道。

4.根据权利要求3所述的稳定域确定方法，其特征在于，所述水轮机的模型表示为：

式中，

m_t为水轮机力矩偏差相对值；q_t为水轮机流量偏差相对值；x为转速偏差相对值，h_t为水轮机水头偏差相对值；y为导叶开度偏差相对值；e_y为水轮机力矩对导叶开度传递系数；e_x为水轮机力矩对转速传递系数；e_h为水轮机力矩对工作水头传递系数；e_qy为水轮机流量对导叶开度传递系数；e_qx为水轮机流量对转速传递系数；e_qh为水轮机流量对工作水头传递系数。

5.根据权利要求4所述的稳定域确定方法，其特征在于，所述发电机组的微分方程表示为：

式中，T_ae为等效机组惯性时间常数，m_te为机组等效主动力矩的偏差相对值，m_g0为等效负载力矩偏差相对值，e_ne为等效发电机负载自调节系数；k_i为第i台机组的额定容量占比，T_ai为第i台机组的惯性时间常数，e_g为发电机负载自调节系数，e_xi为第i台机组力矩对转速传递系数，m_ti为第i台机组的主动力矩偏差相对值，n表示一共有n台机组。

6.根据权利要求5所述的稳定域确定方法，其特征在于，当所述机组为两台时，所述高阶状态空间模型为：

其中，x＝[x,y₁,y′₁,y,y',x_i,x′_i,x_d,x'_d,h_x1,h_x,h'_x1,h'_x,q₂₁,q₂,q'₂₁,q'₂,h_b1,h_b,h′_b1,h′_b,h_p1,h_p,q_c]为状态变量，式中，没有上标的变量为第一台机组单元的变量，有上标“'”的变量表示第二台机组单元的变量；

a_1,1＝-e_ne/T_ae,a_1,4＝ke_y/T_ae,a_1,5＝k'e′_y/T_ae,a_1,11＝-ke_h/T_ae；

a_1,13＝-k'e'_h/T_ae,a_1,19＝ke_h/T_ae,a_1,21＝k'e'_h；

a_2,1＝-K_p/T_y1,a_2,2＝-(1+K_pb_p)/T_y1,a_2,6＝1/T_y1,a_2,8＝1/T_y1；

a_3,1＝-K'_p/T′_y1,a_3,3＝-(1+K'_pb'_p)/T′_y1,a_3,9＝1/T′_y1,a_3,12＝1/T′_y1 a_4,2＝1/T_y,a_4,4＝-1/T_y；

a_5,3＝1/T′_y,a_5,5＝-1/T′_y；

a_6,1＝-K_i,a_6,2＝-K_ib_p；a_7,1＝-K′_i,a_7,3＝-K′_ib'_p；

a_8,1＝(K_de_ne)/(T_dT_ae),a_8,4＝-(kK_de_y)/(T_dT_ae),a_8,5＝-(k'K_de'_y)/(T_dT_ae),a_8,8＝-1/T_d；

a_8,11＝(kK_de_h)/(T_dT_ae),a_8,13＝(k'K_de'_h)/(T_dT_ae),a_8,19＝(kK_de_h)/(T_dT_ae),a_8,21＝-(k'K_de'_h)/(T_dT_ae)；

a_9,1＝(K'_de_ne)/(T′_dT_ae),a_9,4＝-(kK'_de_y)/(T′_dT_ae),a_9,5＝-(k'K'_de'_y)/(T_d'T_ae),a_9,9＝-1/T′_d；

a_9,11＝(kK'_de_h)/(T′_dT_ae),a_9,13＝(k'K'_de'_h)/(T′_dT_ae),a_9,19＝-(kK'_de_h)/(T′_dT_ae),a_9,21＝(k'K'_de'_h)/(T′_dT_ae)；

a_11,1＝8h_w2e_qx/T_r2,a_11,4＝8h_w2e_qy/T_r2,a_11,10＝1,a_11,11＝-8h_w2e_qh/T_r2,a_11,19＝8h_w2e_qh/T_r2；

a_13,1＝8h'_w2e'_qx/T′_r2,a_13,5＝8h'_w2e'_qy/T_r2,a_13,12＝1,a_13,13＝-8h'_w2e'_qh/T′_r2,a_13,21＝8h'_w2e'_qh/T′_r2；

a_15,14＝1,a_15,23＝-2/(h_w2T_r2)；

a_17,16＝1,a_17,23＝-2/(h'_w2T′_r2)；

a_19,1＝-8h_w1e_qx/T_r1,a_19,4＝-8h_w1e_qy/T_r1,a_19,11＝8h_w1e_qh/T_r1,a_19,18＝1,a_19,19＝-8h_w1e_qh/T_r1；

a_21,1＝-8h'_w1e'_qx/T′_r1,a_21,5＝-8h'_w1e'_qy/T′_r1,a_21,13＝8h'_w1e'_qh/T′_r1,a_21,20＝1,a_21,21＝-8h'_w1e'_qh/T′_r1；

a_23,21＝8h_w3/T_r3,a_23,17＝8h_w3/T_r3,a_23,22＝1,a_23,23＝-4a₃/T_r3,a_23,24＝-8h_w3/T_r3；

a_24,15＝8h_w3T_c/(T_r3T_q),a_24,17＝8h_w3T_c/(T_r3T_q),a_24,22＝T_c/T_q；

a_24,23＝-4a₃T_c/(T_r3T_q),a_24,24＝-(8h_w3T_c+T_r3)/(T_r3T_q)；

B为输入矩阵，所述系统矩阵为：

7.根据权利要求1至6任一项所述的稳定域确定方法，其特征在于，通过变步长和二分法确定所述水轮机调节系统的稳定域，具体为：

当微分系数K_d为预设固定值时，将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中：

Step1：采用变步长的方法循环选取比例系数K_p和积分系数K_i的取值为(0，K_{i_0})带入系统矩阵中计算其特征值以判断系统稳定性，进而求得PI临界稳定线的起点(0，K_{i_0})；其中，坐标原点(0，0)到(0，K_{i_0})之间的PI临界稳定线为第一种稳定边界，K_{i_0}为积分系数K_i在第一种稳定边界的最大值；

Step2：按照固定的ΔK_p步长，利用二分法在所述临界稳定线上的上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，具体步骤为：

Step2.1：给定计算精度e,K_p搜索步长dK_p＝0.01，搜索邻域的角度θ范围为初始范围[l,u]初始l＝-pi/2，u＝pi/2；K_i的基准值是K_{i_0}；

Step2.2：在[l，u]区间取中点θ₁；

Step2.3：在θ与dK_p构成三角函数关系：K_i＝dK_p*tanθ₁，计算并判断(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)对应的系统状态；

Step2.4：若系统稳定将θ₁作为θ区间的下限即l＝θ₁，若系统不稳定则将θ₁作为θ区间的上限即u＝θ₁，重复Step2.2～Step2.4，直到|u-l|<e；

Step2.5：将Step2.4满足精度条件的PI参数即(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)作为新Step2.1的基准点，重复Step2.1～Step2.5连续求取新的临界稳定参数，直到K_i<0停止；

Step3：采用变步长的方法循环选取的(K_p，0)带入系统矩阵中计算特征值以判断系统稳定性进而，求得PI临界稳定线的终点(K_{p_end}，0)；其中，坐标原点(0，0)到(K_{p_end}，0)之间的PI临界稳定线为第一种稳定边界，K_{p_end}为比例系数K_p在第一种稳定边界的最大值；

将Step1～Step3求取的PI临界稳定线的数据结合起来，其内部的区域属于所求的稳定域；

将比例系数K_p和积分系数K_i中任一系数设为预设固定值，参见上述Step1～Step3的方法求取微分系数K_d的上限；

改变微分系数K_d的预设固定值，参见上述Step1～Step3求取不同微分系数K_d下的稳定域，将其组合得到系统的稳定域。

8.根据权利要求7所述的稳定域确定方法，其特征在于，所述步骤Step1具体包括如下步骤：

Step1.1：以坐标原点(0，0)为基准点，给定初始步长dk_i＝0.1，在纵轴上找到下一个点判断系统状态，如果稳定则将系统的将步长dk_i＝5*dk_i，并以该点纵坐标为基准更新下一个纵坐标；如果系统不稳定，则将系统步长变为dk_i＝0.2*dk_i，以上一个稳定的点为基础更新纵坐标；

Step1.2：重复Step1.1的更新步骤，直到两次纵坐标的距离小于设定的精度，此时的纵坐标值即为PI临界稳定线的起点(0，K_{i_0})；

所述步骤Step3具体包括如下步骤：

Step3.1：以坐标原点(0，0)为基准点，给定初始步长dk_p＝0.1，在横轴上找到下一个点判断系统状态，如果稳定则将系统的将步长dk_p＝5*dk_p，并以该点横坐标为基准更新下一个横坐标；如果系统不稳定，则将系统步长变为dk_p＝0.2*dk_p，以上一个稳定的点为基础更新横坐标；

Step3.2：重复Step3.1的更新步骤，直到两次横坐标的距离小于设定的精度，此时的横坐标值即为PI临界稳定线的终点(K_{p_end}，0)。

9.一种多机共管道水轮机调节系统稳定域确定系统，其特征在于，包括：

系统模型确定单元，用于确定多机共管道水轮机调节系统高阶状态空间模型；所述水轮机调节系统包括：共管道的多台机组，每台机组包括：PID型调速器、引水系统、水轮机以及发电机四个模块，所述高阶状态空间模型用于通过状态方程表示多机共管道水轮机调节系统的状态变量与各个模块参数组成的系统矩阵之间的关系；

特征值确定单元，用于确定所述系统矩阵的特征值；当所述特征值的实部均小于零时，所述水轮机调节系统处于稳定状态；

稳定域确定单元，用于通过变步长和二分法确定所述水轮机调节系统的稳定域；所述稳定域为能够使所述水轮机调节系统处于稳定状态的所有PID控制参数的集和；所述PID控制参数包括比例系数K_p、积分系数K_i以及微分系数K_d；其中，采用变步长方式确定所述稳定域的第一种稳定边界，所述第一种稳定边界为比例系数K_p、积分系数K_i以及微分系数K_d三个系数中一个系数为预设固定值，另外两个系数一个为0，另一个系数从0到最大值变化；采用二分法确定所述稳定域的第二种稳定边界，所述第二种稳定边界为所述三个系数中一个系数为预设固定值，另外两个系数均从0到最大值变化。

10.根据权利要求9所述的稳定域确定系统，其特征在于，所述稳定域确定单元通过变步长和二分法确定所述水轮机调节系统的稳定域，具体为：

当微分系数K_d为预设固定值时，将K_p-K_i稳定域表示在以K_p为横坐标，K_i为纵坐标的直角坐标系中：

Step1：采用变步长的方法循环选取比例系数K_p和积分系数K_i的取值为(0，K_{i_0})带入系统矩阵中计算其特征值以判断系统稳定性，进而求得PI临界稳定线的起点(0，K_{i_0})；其中，坐标原点(0，0)到(0，K_{i_0})之间的PI临界稳定线为第一种稳定边界，K_{i_0}为积分系数K_i在第一种稳定边界的最大值；

Step2：按照固定的ΔK_p步长，利用二分法在所述临界稳定线上的上个临界稳定点(K_p,K_i)的邻域内求解下一个临界稳定点(K_p+ΔK_p,K_i’)，具体步骤为：

Step2.1：给定计算精度e,K_p搜索步长dK_p＝0.01，搜索邻域的角度θ范围为初始范围[l,u]初始l＝-pi/2，u＝pi/2；K_i的基准值是K_{i_0}；

Step2.2：在[l，u]区间取中点θ₁；

Step2.3：在θ与dK_p构成三角函数关系：K_i＝dK_p*tanθ₁，计算并判断(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)对应的系统状态；

Step2.4：若系统稳定将θ₁作为θ区间的下限即l＝θ₁，若系统不稳定则将θ₁作为θ区间的上限即u＝θ₁，重复Step2.2～Step2.4，直到|u-l|<e；

Step2.5：将Step2.4满足精度条件的PI参数即(K_p＝K_{p_0}+dK_p,K_i＝dK_p*tanθ₁)作为新Step2.1的基准点，重复Step2.1～Step2.5连续求取新的临界稳定参数，直到K_i<0停止；

Step3：采用变步长的方法循环选取的(K_p，0)带入系统矩阵中计算特征值以判断系统稳定性进而，求得PI临界稳定线的终点(K_{p_end}，0)；其中，坐标原点(0，0)到(K_{p_end}，0)之间的PI临界稳定线为第一种稳定边界，K_{p_end}为比例系数K_p在第一种稳定边界的最大值；

将Step1～Step3求取的PI临界稳定线的数据结合起来，其内部的区域属于所求的稳定域；

将比例系数K_p和积分系数K_i中任一系数设为预设固定值，参见上述Step1～Step3的方法求取微分系数K_d的上限；

改变微分系数K_d的预设固定值，参见上述Step1～Step3求取不同微分系数K_d下的稳定域，将其组合得到系统的稳定域。

Images