CN116300426A - 一种带监督机制的推力分配方法 - Google Patents

Abstract

本发明涉及动力定位船舶推力分配技术领域，具体涉及一种带监督机制的推力分配方法，包括：(1)初始化船舶各类参数；(2)判别是否执行最优扇区搜寻，是则运用PSO算法进行最优扇区求解，否则保留前一时刻最优扇区解；(3)判别推进器所处扇区是否与前一时刻选择扇区相一致，是则执行步骤(4)，否则执行步骤(5)；(4)运用非扇区切换的SQP算法进行推力分配求解，然后执行步骤(6)；(5)运用扇区切换的SQP算法进行推力分配求解；(6)运用监督与切换机制进行扇区选择；(7)输出推力分配结果，更新船舶状态参数，跳至步骤(2)进行下一个周期推力分配计算。本发明从全局优化角度进行推力分配求解，保证船舶推力分配精度并降低能耗。

Description

一种带监督机制的推力分配方法

技术领域

本发明涉及动力定位船舶推力分配技术领域，具体涉及一种带监督机制的推力分配方法。

背景技术

船舶动力定位系统作为一种高技术含量的船舶装备已经在各种水上领域特别是海洋领域被广泛应用，诸如熟知的981钻井石油平台、填海“神器”挖泥船等等。推力分配作为船舶动力定位系统的关键环节发挥着至关重要的作用，在控制系统和推进系统之间起到“桥梁”作用，其处理结果直接影响船舶的操纵性、机动性以及能耗等性能。推力分配方法有许多种比如推进器分组法、滤波和奇异值分解法、序列二次规划法(SQP)以及近些年兴起的群智能算法等等。本发明从全局优化角度出发，在综合机械磨损、奇异性和偏差等性能指标前提下降低船舶能耗，提出一种PSO和SQP相融合的带监督机制的推力分配方法。

发明内容

针对上述现有技术中存在的缺点，本发明提供了一种带监督机制的推力分配方法，能够在综合考虑推力分配各项性能指标情况下降低能耗、获取更加可靠及更加优化的推力分配结果。

为了实现上述目的，本发明采用的技术方案如下：

一种带监督机制的推力分配方法，包括如下步骤：

步骤1、初始化船舶各类参数：结构参数、状态参数、物理条件约束参数、权值参数、扇区参数等；

步骤2、判别是否执行最优扇区搜寻，是则运用PSO算法进行最优扇区求解，否则保留前一时刻最优扇区解；

执行最优扇区搜寻条件为：首次执行推力分配或当前需求推力和力矩与上一时刻的推力和力矩不相等；

PSO算法进行最优扇区求解实施方法如下：

步骤2.1、初始化参数，其中粒子选取推力分配方程等式约束项的齐次线性方程组B_eT_e＝0的基础解系数；

T_e为船体坐标系下X方向和Y方向两个相互正交的力，(l_xn,l_yn)为第n个推进器在船体坐标系中的位置；

步骤2.2、以能耗项和偏差项为代价函数计算粒子适应度值，并计算全局最优值，该步骤中推进器角度在禁区范围内处理方法为：计算该角度与该禁区极值的距离，靠近禁区极大值则取该值为推进器角度值，否则取禁区极小值为推进器角度值；

步骤2.3、运用带权重系数的PSO算法进行种群进化，粒子更新；

步骤2.4、判别迭代次数是否达设定值，是执行步骤2.5，否执行步骤2.3；

步骤2.5、对所求最优推力分配解的推进器角度进行扇区判别，具体实施方法如下：首先判别角度值是否处于前一时刻选择扇区内，是则判定最优扇区取前一时刻选择扇区，否则判别角度值是否处于前一时刻之前选择扇区内，是则判定最优扇区取前一时刻之前选择扇区，否则按照扇区序数从小到大逐一判别进行最优扇区取值；

步骤3、判别推进器所处扇区是否与前一时刻选择扇区相一致，是则执行步骤4，否则执行步骤5；

步骤4、运用非扇区切换的SQP算法进行推力分配求解，然后执行步骤6，非扇区切换的SQP算法具体实施方法如下：

构建二次型推力分配数学模型为：

式中H是包含能耗项(采用一阶泰勒级数展开)、偏差项、机械磨损项和奇异项(采用二阶泰勒级数展开)的二次项权系数矩阵，f是对应的一次项权系数向量，X是包含推力变化量、方位角变化量和松弛变量的优化变量，A_eq和b_eq分别为等式约束矩阵和向量，L_b和U_b为包含推力变化量、方位角变化量和松弛变量的优化变量边界；

b_eq＝τ-B(α_P)T_P

ΔST_imin＝max(ΔT_imin,T_imin-T_iP)，ΔST_imax＝min(ΔT_imax,T_imax-T_iP)

ΔSα_imax＝min(Δα_imax,α_imax(S_now)-α_iP)，ΔSα_imin＝max(Δα_imin,α_imin(S_now)-α_iP)；

其中B(α_P)是前一时刻推进器配置控制矩阵，α_P和T_P分别为前一时刻方位角和推力，E₃为三阶单位矩阵，τ为控制系统输出的推力和力矩指令，ΔST_imin和ΔT_imin分别是第i个推进器当前和自身物理条件允许的推力最小变化量，ΔST_imax和ΔT_imax分别是第i个推进器当前和自身物理条件允许的推力最大变化量，ΔSα_imin和ΔSα_imax分别是第i个推进器允许的方位角最小和最大变化量，Δα_imax和Δα_imin分别是第i个推进器自身物理条件允许的方位角最大和最小变化量，α_imax(S_now)和α_imin(S_now)为第i个推进器前一时刻选择扇区的边界，α_iP和T_iP分别为第i个推进器当前角度和推力，S_now为前一时刻选择扇区值；

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤5、运用扇区切换的SQP算法进行推力分配求解，具体实施方法如下：

构建二次型推力分配数学模型为：

式中H_S是包含能耗项(采用一阶泰勒级数展开)和偏差项的二次项权系数矩阵，f_S是对应的一次项权系数向量，X_S是包含推力变化量和松弛变量的优化变量，A_Seq和b_Seq分别为等式约束矩阵和向量，L_Sb和U_Sb为包含推力变化量和松弛变量的优化变量边界，其中A_Seq和b_Seq计算方法如下：

A_Seq＝[B(α_P+Δα)| E₃]，b_Seq＝τ-B(α_P+Δα)T_P

其中B(α_P+Δα)为此刻的推进器配置控制矩阵，Δα为方位角变化量；

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤6、运用监督与切换机制进行扇区选择，该机制实施思路为：前一时刻选择扇区每次更新只有两种可能，一是保持不变，二是被最优扇区取代，而要被最优扇区取代只要满足如下三个条件中一个：(1)最优扇区合力和力矩偏差比当前扇区小于一定值，且保持一段时间；(2)最优扇区功率比当前扇区小于一定值，最优扇区合力和力矩偏差比当前扇区不大于一定值，且保持一段时间；(3)最优扇区与前一时刻之前选择扇区相等且切换时间不超过一定值；

步骤7、输出推力分配结果，更新船舶状态参数，跳至步骤2进行下一个周期推力分配计算。

优选的，所述步骤2中，推进器角度在禁区范围内处理具体计算方法如下：

设推进器在禁区内的角度为α，禁区范围为[α_strmin,α_strmax]，则α处理结果为：

优选的，所述步骤5中，SQP算法进行推力分配求解中方位角变化量具体计算方法为：

优选的，所述步骤6中，运用监督与切换机制进行扇区选择具体实施计算方法如下：

假设前一时刻之前选择的扇区为S_ago，当前最优扇区为S_best，当前选择的扇区为S_choose，则：

t是指S_best和S_now连续不等累加的时间，t_lim是指停留时间切换的最小切换时间，t₀是为了防止瞬态干扰导致的切换所设置的最大返回切换时间，P_best是全局最优消耗功率值，P_SQP是扇区切换或非扇区切换的SQP算法求解所得消耗功率值，e_SQP是扇区切换或非扇区切换的SQP算法求解结果偏差的平方和，e_best是全局最优求解结果偏差的平方和，e_S为偏差设定值，P₁和P₂为设定的百分数。

与现有技术相比，本发明具有以下有益效果：

1、本发明首先利用PSO算法的全局寻优优势进行最优扇区求解，然后根据扇区逻辑选择扇区切换的SQP算法或非扇区切换的SQP算法进行推力分配求解，最后根据监督与切换机制进行扇区选择为下次推力分配求解指明方向，实验表明该算法能满足工程应用实时性需求又能通过扇区切换灵活优化推力分配各项指标。

2、本发明所提出的推力分配方法能够在考虑能耗和偏差指标情况下获取全局最优解及最优扇区组合，在此基础上根据监督与切换机制进行扇区选择，并能根据扇区逻辑合理选择SQP推力分配算法，进而获取更加可靠及更加优化的推力分配结果。

附图说明

图1是本发明的流程示意图；

图2是本发明纵向力输出曲线图；

图3是本发明横向力输出曲线图；

图4是本发明力矩输出曲线图；

图5是本发明推力输出曲线图；

图6是本发明方位角输出曲线图；

图7是本发明推进器消耗功率输出曲线图。

具体实施方式

下面结合附图1-7和船模实例对本发明的技术方案进行详细说明：

一种带监督机制的推力分配方法，包括如下步骤：

步骤1、初始化船舶各类参数：结构参数、状态参数、物理条件约束参数、权值参数、扇区参数等，具体参数如下所示：

表1推进器的部分参数：

推进器	位置(m)	Tmax/Tmin(N)	ΔTmax/ΔTmin(N)	Δαmax/Δαmin(rad)	Kpt
						1号	(-0.875,-0.123)	21.9/0	4.38/-4.38	0.17/-0.17	1.5844
2号	(-0.875,0.123)	21.9/0	4.38/-4.38	0.17/-0.17	1.5844
						3号	(0.55,0)	10/-10	2/-2	0/0	41.919

S_ago，S_best，S_choose，S_now的初始值皆为1，初始推力大小皆为0，推进器1～3的初始方位角分别为0.01、-0.01和1.57弧度，t_lim＝4，t₀＝2，推进器1的扇区分为：[-1.22,1.9]，[-0.349,2.775]，[-2.618,-1.92]，推进器2的扇区分为：[-1.9,1.22]，[-2.775,0.349]，[1.92,2.618]，推进器3始终为1.57弧度，故共有9种扇区组合；

步骤2、判别是否执行最优扇区搜寻，是则运用PSO算法进行最优扇区求解，否则保留前一时刻最优扇区解，执行最优扇区搜寻条件为：首次执行推力分配或当前需求推力和力矩与上一时刻的推力和力矩不相等，PSO算法进行最优扇区求解实施方法如下：

步骤2.1、初始化参数；

步骤2.2、以能耗项和偏差项为代价函数计算粒子适应度值，并计算全局最优值，该步骤中推进器角度在禁区范围内处理方法为：计算该角度与该禁区极值的距离，靠近禁区极大值则取该值为推进器角度值，否则取禁区极小值为推进器角度值；

步骤2.3、运用带权重系数的PSO算法进行种群进化，粒子更新；

步骤2.4、判别迭代次数是否达设定值，是执行步骤2.5，否执行步骤2.3；

步骤2.5、对所求最优推力分配解的推进器角度进行扇区判别，具体实施方法如下：首先判别角度值是否处于前一时刻选择扇区内，是则判定最优扇区取前一时刻选择扇区，否则判别角度值是否处于前一时刻之前选择扇区内，是则判定最优扇区取前一时刻之前选择扇区，否则按照扇区序数从小到大逐一判别进行最优扇区取值；

其中，推进器角度在禁区范围内处理具体计算方法如下：

设推进器在禁区内的角度为α，禁区范围为[α_strmin,α_strmax]，则α处理结果为：

步骤3、判别推进器所处扇区是否与前一时刻选择扇区相一致，是则执行步骤4，否则执行步骤5；

步骤4、运用非扇区切换的SQP算法进行推力分配求解，具体实施方法如下：

构建二次型推力分配数学模型为：

式中H是包含能耗项(采用一阶泰勒级数展开)、偏差项、机械磨损项和奇异项(采用二阶泰勒级数展开)的二次项权系数矩阵，f是对应的一次项权系数向量，X是包含推力变化量、方位角变化量和松弛变量的优化变量，A_eq和b_eq分别为等式约束矩阵和向量，L_b和U_b为包含推力变化量、方位角变化量和松弛变量的优化变量边界；

b_eq＝τ-B(α_P)T_P

ΔST_imin＝max(ΔT_imin,T_imin-T_iP)，ΔST_imax＝min(ΔT_imax,T_imax-T_iP)

ΔSα_imax＝min(0.17,α_imax(S_now)-α_iP)，ΔSα_imin＝max(-0.17,α_imin(S_now)-α_iP)；

其中B(α_P)是前一时刻推进器配置控制矩阵，α_P和T_P分别为前一时刻方位角和推力，E₃为三阶单位矩阵，τ为控制系统输出的推力和力矩指令，ΔST_imin和ΔT_imin分别是第i个推进器当前和自身物理条件允许的推力最小变化量，ΔST_imax和ΔT_imax分别是第i个推进器当前和自身物理条件允许的推力最大变化量，ΔSα_imin和ΔSα_imax分别是第i个推进器允许的方位角最小和最大变化量，α_imax(S_now)和α_imin(S_now)为第i个推进器前一时刻选择扇区的边界，α_iP和T_iP分别为第i个推进器当前角度和推力，S_now为前一时刻选择扇区值；

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤5、运用扇区切换的SQP算法进行推力分配求解，具体实施方法如下：

构建二次型推力分配数学模型为：

式中H_S是包含能耗项(采用一阶泰勒级数展开)和偏差项的二次项权系数矩阵，f_S是对应的一次项权系数向量，X_S是包含推力变化量和松弛变量的优化变量，A_Seq和b_Seq分别为等式约束矩阵和向量，L_Sb和U_Sb为包含推力变化量和松弛变量的优化变量边界，其中A_Seq和b_Seq计算方法如下：

A_Seq＝[B(α_P+Δα)| E₃]，b_Seq＝τ-B(α_P+Δα)T_P

b_Seq＝τ-B(α_P+Δα)T_P

其中B(α_P+Δα)为此刻的推进器配置控制矩阵，Δα为方位角变化量其具体计算方法如下：

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤6、运用监督与切换机制进行扇区选择，具体实施方法如下：

假设前一时刻之前选择的扇区为S_ago，当前最优扇区为S_best，当前选择的扇区为S_choose，则：

t是指S_best和S_now连续不等累加的时间，P_best是全局最优消耗功率值，P_SQP是扇区切换或非扇区切换的SQP算法求解所得消耗功率值，e_SQP是扇区切换或非扇区切换的SQP算法求解结果偏差的平方和，e_best是全局最优求解结果偏差的平方和；

步骤7、输出推力分配结果，更新船舶状态参数，跳至步骤2进行下一个周期推力分配计算。

结果说明：

图2至图4中曲线1、2分别对应设定值和实际输出值，图5至图6中曲线1、2、3分别对应推进器1、2、3的相应输出值，从图2、图3和图4可以看出，采取本发明推力分配方法输出的实际纵向力、横向力和力矩与控制系统期望的纵向力、横向力和力矩几乎相等，达到了非常好的推力分配效果，推力在采样周期101S时刻产生突变(该突变使得最优解方向正好与原最优方位相反，最优扇区发生改变需要推进器大范围旋转)情况下，推进器能以最快速度达到新的稳定态，在过渡过程中一部分时间段内推力解是不能满足期望值的，为了降低能耗各个推进器的推力降低到0附近，这也说明本算法在不同情况下都能尽量把能耗降到最低，图6的角度变化趋势可以看出推力分配结果避免了奇异性，同时两个推进器角度几乎相反方向等值变化有利于能耗的最小化，图7可以看出除了在过渡过程中功率方式短暂的突变，其它时间段内变化平缓且和理论最优值非常靠近(不考虑奇异性时两种情况下最优功率分别为8.96W和3.1688W，而实际输出为9.048W和3.512W)，这也说明本发明的推力分配方法可以降低功耗，有利于节能减排。

综上所述，本发明所提出的推力分配方法能够在综合考虑推力分配各项性能指标情况下灵活的优化各种性能指标、降低能耗进而获取更加可靠及更加优化的推力分配结果。

以上对本发明的一个实施例进行了详细说明，但所述内容仅为本发明的较佳实施例，不能被认为用于限定本发明的实施范围。凡依本发明申请范围所作的均等变化与改进等，均应仍归属于本发明的专利涵盖范围之内。

Claims (4)

1.一种带监督机制的推力分配方法，其特征在于：包括如下步骤：

步骤1、初始化船舶各类参数：结构参数、状态参数、物理条件约束参数、权值参数、扇区参数；

步骤2、判别是否执行最优扇区搜寻，是则运用PSO算法进行最优扇区求解，否则保留前一时刻最优扇区解；

执行最优扇区搜寻条件为：首次执行推力分配或当前需求推力和力矩与上一时刻的推力和力矩不相等；

PSO算法进行最优扇区求解实施方法如下：

步骤2.1、初始化参数，其中粒子选取推力分配方程等式约束项的齐次线性方程组B_eT_e＝0的基础解系数；

T_e为船体坐标系下X方向和Y方向两个相互正交的力，(l_xn,l_yn)为第n个推进器在船体坐标系中的位置；

步骤2.2、以能耗项和偏差项为代价函数计算粒子适应度值，并计算全局最优值，该步骤中推进器角度在禁区范围内处理方法为：计算该角度与该禁区极值的距离，靠近禁区极大值则取该值为推进器角度值，否则取禁区极小值为推进器角度值；

步骤2.3、运用带权重系数的PSO算法进行种群进化，粒子更新；

步骤2.4、判别迭代次数是否达设定值，是执行步骤2.5，否执行步骤2.3；

步骤2.5、对所求最优推力分配解的推进器角度进行扇区判别，具体实施方法如下：首先判别角度值是否处于前一时刻选择扇区内，是则判定最优扇区取前一时刻选择扇区，否则判别角度值是否处于前一时刻之前选择扇区内，是则判定最优扇区取前一时刻之前选择扇区，否则按照扇区序数从小到大逐一判别进行最优扇区取值；

步骤3、判别推进器所处扇区是否与前一时刻选择扇区相一致，是则执行步骤4，否则执行步骤5；

步骤4、运用非扇区切换的SQP算法进行推力分配求解，然后执行步骤6，非扇区切换的SQP算法具体实施方法如下：

构建二次型推力分配数学模型为：

式中H是包含能耗项、偏差项、机械磨损项和奇异项的二次项权系数矩阵，f是对应的一次项权系数向量，X是包含推力变化量、方位角变化量和松弛变量的优化变量，A_eq和b_eq分别为等式约束矩阵和向量，L_b和U_b为包含推力变化量、方位角变化量和松弛变量的优化变量边界；

b_eq＝τ-B(α_P)T_P

ΔST_imin＝max(ΔT_imin,T_imin-T_iP)，ΔST_imax＝min(ΔT_imax,T_imax-T_iP)

ΔSα_imax＝min(Δα_imax,α_imax(S_now)-α_iP)，ΔSα_imin＝max(Δα_imin,α_imin(S_now)-α_iP)

其中B(α_P)是前一时刻推进器配置控制矩阵，α_P和T_P分别为前一时刻方位角和推力，E₃为三阶单位矩阵，τ为控制系统输出的推力和力矩指令，ΔST_imin和ΔT_imin分别是第i个推进器当前和自身物理条件允许的推力最小变化量，ΔST_imax和ΔT_imax分别是第i个推进器当前和自身物理条件允许的推力最大变化量，ΔSα_imin和ΔSα_imax分别是第i个推进器允许的方位角最小和最大变化量，Δα_imax和Δα_imin分别是第i个推进器自身物理条件允许的方位角最大和最小变化量，α_imax(S_now)和α_imin(S_now)为第i个推进器前一时刻选择扇区的边界，α_iP和T_iP分别为第i个推进器当前角度和推力，S_now为前一时刻选择扇区值；

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤5、运用扇区切换的SQP算法进行推力分配求解，具体实施方法如下：构建二次型推力分配数学模型为：

式中H_S是包含能耗项和偏差项的二次项权系数矩阵，f_S是对应的一次项权系数向量，X_S是包含推力变化量和松弛变量的优化变量，A_Seq和b_Seq分别为等式约束矩阵和向量，L_Sb和U_Sb为包含推力变化量和松弛变量的优化变量边界，其中A_Seq和b_Seq计算方法如下：

A_Seq＝[B(α_P+Δα)|E₃]，b_Seq＝τ-B(α_P+Δα)T_P

其中B(α_P+Δα)为此刻的推进器配置控制矩阵，Δα为方位角变化量；

根据构建的二次型推力分配数学模型利用SQP算法进行推力求解；

步骤6、运用监督与切换机制进行扇区选择，该机制实施思路为：前一时刻选择扇区每次更新只有两种可能，一是保持不变，二是被最优扇区取代，而要被最优扇区取代只要满足如下三个条件中一个：(1)最优扇区合力和力矩偏差比当前扇区小于一定值，且保持一段时间；(2)最优扇区功率比当前扇区小于一定值，最优扇区合力和力矩偏差比当前扇区不大于一定值，且保持一段时间；(3)最优扇区与前一时刻之前选择扇区相等且切换时间不超过一定值；

步骤7、输出推力分配结果，更新船舶状态参数，跳至步骤2进行下一个周期推力分配计算。

2.根据权利要求1所述的一种带监督机制的推力分配方法，其特征在于：所述步骤2中，推进器角度在禁区范围内处理具体计算方法如下：

设推进器在禁区内的角度为α，禁区范围为[α_strmin,α_strmax]，则α处理结果为：

3.根据权利要求1所述的一种带监督机制的推力分配方法，其特征在于：所述步骤5中，SQP算法进行推力分配求解中方位角变化量具体计算方法为：

4.根据权利要求1所述的一种带监督机制的推力分配方法，其特征在于：所述步骤6中，运用监督与切换机制进行扇区选择具体实施计算方法如下：

假设前一时刻之前选择的扇区为S_ago，当前最优扇区为S_best，当前选择的扇区为S_choose，则：

t是指S_best和S_now连续不等累加的时间，t_lim是指停留时间切换的最小切换时间，t₀是为了防止瞬态干扰导致的切换所设置的最大返回切换时间，P_best是全局最优消耗功率值，P_SQP是扇区切换或非扇区切换的SQP算法求解所得消耗功率值，e_SQP是扇区切换或非扇区切换的SQP算法求解结果偏差的平方和，e_best是全局最优求解结果偏差的平方和，e_S为偏差设定值，P₁和P₂为设定的百分数。

Images