一种动力定位推力系统智能推力分配方法.pdf

本发明涉及一种动力定位推力系统智能推力分配方法，包括以下步骤：(1)推力系统配置参数输入；(2)数据准备与存储；(3)系统对准与优化建模；(4)优化问题求解；(5)推力分配决策输出。与现有技术相比，本发明具有实时性好、人机交互友好、人员工作量小、自动化程度高、决策效果好、能实现较好的推力分配等优点。

权利要求书
1.  一种动力定位推力系统智能推力分配方法，其特征在于，包括以下步骤：
(1)推力系统配置参数输入；
(2)数据准备与存储；
(3)系统对准与优化建模；
(4)优化问题求解；
(5)推力分配决策输出。

2.  根据权利要求1所述的智能推力分配方法，其特征在于，所述的推力系统配置参数输入具体为：
11)输入推力系统中推力设备的个数p；
12)对p个推力设备通过窗口方式进行参数设置，参数设置具体包括：
①是否可变方向；
②设备安装位置，即船体坐标系下的坐标值，存储为变量lk，x，lk，y；
③产生推力的上限值，存储为变量Tkmax；
④安装角度，即船体坐标系下的角度值，对于①中参数为“否”的推力设备，需要输入推力设备产生正推力方向的推力角度值αk；对于①中参数为“是”的推力设备，需要输入推力设备产生正推力方向的角度范围，存储为变量αkmin和αkmax；
⑤选择推力设备灵敏程度，包括高档、中档、低档和自定义，存储为变量εk；
⑥优化参数选择，包括默认和自定义，存储为变量Hk，Mk，Qk。

3.  根据权利要求2所述的智能推力分配方法，其特征在于，所述的数据准备与存储具体为：
21)自行根据参数将①中参数为“是”的推力设备进行编号为1，2，…，pr，将①中参数为“否”的推力设备进行编号为pr+1，…，p，同时将各推力设备变量与编号对应；
22)计算并存储推力系数矩阵B：B=[Br，Bf]
其中：
Br=10...1001...01-l1,yl1,x...-lpr,ylpr,xBf=cosαpr+1...cosαpsinαpr+1...sinαplpr+1...lp]]>
其中，

23)计算并存储4k组变量，mki，nki，其中k＝1，2，…，pr，i=1，2，3，4：
Dk={α|αkmin≤α≤αkmax}，k＝1，2，…，pr，
Ft={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>
Dki=Dk∩Fi，
mki=inf{tanα|α∈Dki}，nki=sup{tanα|α∈Dki}；
24)计算并存储组二次规划变量，每组编号为t，其中将t进行4进制转换：
t=(dprdpr-1…d1)4
将4进制的第(j-1)位存储为dj，其中j=1，2，…，pr；
第t组需计算和存储的约束变量为：

其中：
At=[A1，A2，t，A3，t]T，Bt=[b1，0]T，Ct=[0，C3]T
A1=[E，-E]T，b1=[b11，b12]T
b11=[22T1,...,22Tpr,T(pr+1)max,...,Tpmax]T,b12=-b11;]]>

A2,t=A2,t,1...A2,t,prA3,t=A3,t,1...A3,t,prC3=C3,1...C3,pr]]>
A2temp=(mkie2k-1-e2k)T(-nkie2k-1+e2k)TA2,t,k=A2temp,dk=0,3-A2temp,dk=1,2C3,k=-e2kTe2kT]]>
A3temp1=(-ϵke2k-1-e2k)T(-ϵke2k-1+e2k)TA3temp2=(ϵke2k-1-e2k)T(ϵke2k-1+e2k)TA3,t,k=A3temp1,dk=0,3A3temp2,dk=1,2]]>
其中ej表示单位矩阵E的第j列，其中At、Bt、Ct为计算过程的中间变量；
25)计算并存储优化参数
若优化参数选项为“默认”，则选择默认的Hdefault，Mdefault，Qdefault数值赋值给Hk，Mk，Qk；若优化参数选项为“自定义”，则按照推力设备的输入信息记录Hk，Mk，Qk。
然后存储如下的优化参数矩阵：
H=diag[H1,H1,...,Hpr,Hpr,Hpr+1,Hpr+2,...,Hp]]]>
M=diag[M1,M1,...,Mpr,Mpr,Mpr+1,Mpr+2,...,Mp]]]>
Q=diag[Q1,Q1,...,Qpr,Qpr,Qpr+1,Qpr+2,...,Qp]]]>
。

4.  根据权利要求3所述的智能推力分配方法，其特征在于，所述的系统对准与优化建模具体为：
31)获取推力设备的状态数据，包括每个推力设备当前的输出推力TRk(k=1，2，…，p)以及当前时刻推力输出的方向αRk(k=1，2，…，pr)；
32)计算u0并存储，在系统运行的过程中需要不断刷新：
u0=[TR0cosαR0,TR0sinαR0,...,TRprcosαRpr,TRprsinαRpr,TR(pr+1),TR(pr+2),...,TRp]]]>
33)输入需要的总推力数据τc，其中τc=[Xc，Yc，Nc]T
其中Xc，Yc分别为推力系统在进退、横移方向所需推进力，Nc为艏摇方向所需推进力矩；
34)利用给定数据生成个子可行域上的凸二次规划问题，第t个二次规划问题的标准形式为：
minJt(u，s)=uTHu+(u-u0)TM(u-u0)+sTQs
s.t.    τc=Bu+s
Atu≤Bt+Ctu0
其中为松弛变量，表示控制输出和给定值之间的差值。

5.  根据权利要求4所述的智能推力分配方法，其特征在于，所述的子可行域上的凸二次规划的角度约束范围为0～90°，90°～180°，180°～270°，270°～360°的其中之一。

6.  根据权利要求4所述的智能推力分配方法，其特征在于，所述的优化问题求解具体为：
根据参数生成的个二次规划问题，每个问题会产生最优解ut，st，以及相应的指标值Jt；
比较这个二次规划的结果，选取：
tm=argmint∈IJt,I={1,2,...,4pr}]]>
并将作为最终的决策向量存储。

7.  根据权利要求6所述的智能推力分配方法，其特征在于，所述的推力分配决策输出具体为：
根据决策向量计算每个推进设备的角度和推力数值：
Tk=utmT(E2k-1+E2k)utm,k=1,2,...,prepr+kutm,k=pr+1,...,n]]>
其中ej表示单位矩阵E的第j列；Ej为正交投影矩阵，其中单位矩阵
推进设备角度计算，tm转换为4进制表示：
tm=(dprdpr-1…d1)4
根据条件：
tanαk=YkXk=e2kTutme2k-1Tutmαk∈Fi=dk+1]]>
其中Fi={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>求解出αk的数值；
最终给出的动力分配方案为：
对于可旋转的推进器的动力分配决策为(Tk，αk)，k=1，2，…，pr；对于固定的推进器的动力分配决策为Tk，k=pr+1，…，n。

8.  根据权利要求7所述的智能推力分配方法，其特征在于，所述的动力分配决策输出可选择多种总线输出方式，包括RS485\CAN\PROFIBUS总线模式，使其适应不同类型的动力定位推力系统总线类型。

说明书一种动力定位推力系统智能推力分配方法
技术领域
本发明涉及一种推力系统的推力分配方法，尤其是涉及一种动力定位推力系统智能推力分配方法。
背景技术
动力定位系统中，推力分配是指针对定位系统中的多个推进器，在允许的时间里，同时在考虑诸多物理约束的条件下，计算出安装在平台上的每个推进器应该在什么方向上产生多大的推力，而且找出可能的多个推力和方向组合中最好的那一种分配，并产生相应的执行指令。
目前国内动力定位中推力优化分配方法研究与系统设计的局限性主要体现在如下两个方面：
其一，体现在设计和开发以及维护成本方面。目前国外推力分配开发模式一般是，将推力分配作为动力定位系统的一部分交由企业或研究机构完成，而国内则一般引进国外技术。这种模式的缺点在于：定制周期长，定制成本高，只能适用与具体动力定位对象，设备维护后期系统维护困难，用户很难根据自己的需求进行推力分配等。
其二，体现在推力分配的优化模型建立过程中，由于推进器的约束和灵敏度的限制，传统的优化建模方法使得推力优化分配区域变成了非凸集(T.A.Johansen，T.I.Fossen，S.P.Berge，Constrained nonlinear control allocation with singularity avoidance using sequential quadratic programming，IEEE Transaction on Control Systems Technology，2004，12(1)：211-216.)，采用各种优化算法效果不甚理想，全局最优解的求解困难。
发明内容
本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种实时性好、 人机交互友好、人员工作量小、自动化程度高、决策效果好、能实现较好的推力分配的动力定位推力系统智能推力分配方法。
本发明的目的可以通过以下技术方案来实现：
一种动力定位推力系统智能推力分配方法，其特征在于，包括以下步骤：
(1)推力系统配置参数输入；
(2)数据准备与存储；
(3)系统对准与优化建模；
(4)优化问题求解；
(5)推力分配决策输出。
所述的推力系统配置参数输入具体为：
11)输入推力系统中推力设备的个数p；
12)对p个推力设备通过窗口方式进行参数设置，参数设置具体包括：
①是否可变方向；
②设备安装位置，即船体坐标系下的坐标值，存储为变量lk，x，lk，y；
③产生推力的上限值，存储为变量Tkmax；
④安装角度，即船体坐标系下的角度值，对于①中参数为“否”的推力设备，需要输入推力设备产生正推力方向的推力角度值αk；对于①中参数为“是”的推力设备，需要输入推力设备产生正推力方向的角度范围，存储为变量αkmin和αkmax；
⑤选择推力设备灵敏程度，包括高档、中档、低档和自定义，存储为变量εk；
⑥优化参数选择，包括默认和自定义，存储为变量Hk，Mk，Qk。
所述的数据准备与存储具体为：
21)自行根据参数将①中参数为“是”的推力设备进行编号为1，2，…，pr，将①中参数为“否”的推力设备进行编号为pr+1，…，p，同时将各推力设备变量与编号对应；
22)计算并存储推力系数矩阵B：B=[Br，Bf]
其中：
Br=10...1001...01-l1,yl1,x...-lpr,ylpr,xBf=cosαpr+1...cosαpsinαpr+1...sinαplpr+1...lp]]>
其中，

23)计算并存储4k组变量，mki，nki，其中k＝1，2，…，pr，i＝1，2，3，4：
Dk={α|αkmin≤α≤αkmax}，k=1，2，…，pr，
Ft={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>
Dki=Dk∩Fi，
mki＝inf{tanα|α∈Dki}，nki=sup{tanα|α∈Dki}；
24)计算并存储组二次规划变量，每组编号为t，其中将t进行4进制转换：
t=(dprdpr-1…d1)4
将4进制的第(j-1)位存储为dj，其中j＝1，2，…，pr；
第t组需计算和存储的约束变量为：

其中：
At=[A1，A2，t，A3，t]T，Bt=[b1，0]T，Ct=[0，C3]T
A1=[E，-E]T，b1＝[b11，b12]T
b11=[22T1,···,22Tpr,T(pr+1)max,···,Tpmax]T,b12=-b11;]]>

A2,t=A2,t,1...A2,t,prA3,t=A3,t,1...A3,t,prC3=C3,1...C3,pr]]>
A2temp=(mkie2k-1-e2k)T(-nkie2k-1+e2k)TA2,t,k=A2temp,dk=0,3-A2temp,dk=1,2C3,k=-e2kTe2kT]]>
A3temp1=(-ϵke2k-1-e2k)T(-ϵke2k-1+e2k)TA3temp2=(ϵke2k-1-e2k)T(ϵke2k-1+e2k)TA3,t,k=A3temp1,dk=0,3A3temp2,dk=1,2]]>
其中ej表示单位矩阵E的第j列，其中At、Bt、Ct为计算过程的中间变量；
25)计算并存储优化参数
若优化参数选项为“默认”，则选择默认的Hdefault，Mdefault，Qdefault数值赋值给Hk，Mk，Qk；若优化参数选项为“自定义”，则按照推力设备的输入信息记录Hk，Mk，Qk。
然后存储如下的优化参数矩阵：
H=diag[H1,H1,...,Hpr,Hpr,Hpr+1,Hpr+2,...,Hp]]]>
M=diag[M1,M1,...,Mpr,Mpr,Mpr+1,Mpr+2,...,Mp]]]>
Q=diag[Q1,Q1,...,Qpr,Qpr,Qpr+1,Qpr+2,...,Qp]]]>
。
所述的系统对准与优化建模具体为：
31)获取推力设备的状态数据，包括每个推力设备当前的输出推力TRk(k=1，2，…，p)以及当前时刻推力输出的方向αRk(k=1，2，…，pr)；
32)计算u0并存储，在系统运行的过程中需要不断刷新：
u0=[TR0cosαR0,TR0sinαR0,...,TRprcosαRpr,TRprsinαRpr,TR(pr+1),TR(pr+2),...,TRp]]]>
33)输入需要的总推力数据τc，其中τc=[Xc，Yc，Nc]T
其中Xc，Yc分别为推力系统在进退、横移方向所需推进力，Nc为艏摇方向所需推进力矩；
34)利用给定数据生成个子可行域上的凸二次规划问题，第t个二次规划问题的标准形式为：
minJt(u，s)＝uTHu+(u-u0)TM(u-u0)+sTQs
s.t.    τc=Bu+s
Atu≤Bt+Ctu0
其中为松弛变量，表示控制输出和给定值之间的差值。
所述的子可行域上的凸二次规划的角度约束范围为0～90°，90°～180°，180°～270°，270°～360°的其中之一。
所述的优化问题求解具体为：
根据参数生成的个二次规划问题，每个问题会产生最优解ut，st，以及相应的指标值Jt；
比较这个二次规划的结果，选取：
tm=argmint∈IJt,I={1,2,...,4pr}]]>
并将作为最终的决策向量存储。
所述的推力分配决策输出具体为：
根据决策向量计算每个推进设备的角度和推力数值：
Tk=utmT(E2k-1+E2k)utm,k=1,2,...,prepr+kutm,k=pr+1,...,n]]>
其中ej表示单位矩阵E的第j列；Ej为正交投影矩阵，其中单位矩阵
推进设备角度计算，tm转换为4进制表示：
tm=(dprdpr-1…d1)4
根据条件：
tanαk=YkXk=e2kTutme2k-1Tutmαk∈Fi=dk+1]]>
其中Fi={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>求解出αk的数值；
最终给出的动力分配方案为：
对于可旋转的推进器的动力分配决策为(Tk，αk)，k=1，2，…，pr；对于固定的推进器的动力分配决策为Tk，k=pr+1，…，n。
将Tk，αk存入指定内存区域中。Tk，αk指定内存区域中的数据随方法运行不断刷新，当在总线(CAN\RS485\PROFIBUS)中工作过程中遇到数据请求时，就将内存中Tk，αk的数据打包成CAN\RS485\PROFIBUS数据帧格式按照相应总线的通信协议发送上层或下层指定控制设备。
所述的动力分配决策输出可选择多种总线输出方式，包括RS485\CAN\PROFIBUS总线模式，使其适应不同类型的动力定位推力系统总线类型。
与现有技术相比，本发明具有以下优点：
1、使用本发明方法进行推力分配设计时，用户只需要了解具体动力定位推力系统配置参数，并将参数输入系统，推力分配部分具体实现由搭载算法的计算机或控制系统完成，使得用户不必具备推力分配设计经验和资历，缩短了开发周期和相应成本；同时，对于高要求的用户，可以通过理论分析或测试的方式修改或自定义优化参数，实现推力分配方案的定制，另外，这样的设计，同时也能解决在推力系统某推进装置更换后，与原推力分配设置不兼容的情况，提高了推力分配的灵活性，缩短了维护周期和相应成本。
2、本推力分配方法在解决这个问题的过程中，将原本的非凸集推力优化分配区域进行了分解，转化为有限个凸多面集上，从而将原二次规划问题转化为有限个在子可行域上的凸二次规划问题，再通过这些凸二次规划问题的求解与比较，计算 出全局最优推力分配方案，这种转化对求取推力分配的最优决策是十分必要的，解决了全局最优解求解困难的问题。
附图说明
图1为本发明的智能推力分配方法的流程图；
图2为本发明优化求解流程图；
图3是本发明配置参数输入界面；
图4是本发明针对具体实施案例用户推力系统配置参数输入图。
具体实施方式
下面结合附图和具体实施例对本发明进行详细说明。
实施例
如图1和图2所示，一种动力定位推力系统智能推力分配方法，其特征在于，包括以下步骤：
(1)推力系统配置参数输入；
(2)数据准备与存储；
(3)系统对准与优化建模；
(4)优化问题求解；
(5)推力分配决策输出。
所述的推力系统配置参数输入具体为：
11)输入推力系统中推力设备的个数p；
12)对p个推力设备通过窗口方式进行参数设置，参数设置具体包括：
①是否可变方向；
②设备安装位置，即船体坐标系下的坐标值，存储为变量lk，x，lk，y；
③产生推力的上限值，存储为变量Tkmax；
④安装角度，即船体坐标系下的角度值，对于①中参数为“否”的推力设备，需要输入推力设备产生正推力方向的推力角度值αk；对于①中参数为“是”的推力设备，需要输入推力设备产生正推力方向的角度范围，存储为变量αkmin和αkmax；
⑤选择推力设备灵敏程度，包括高档、中档、低档和自定义，存储为变量εk；
⑥优化参数选择，包括默认和自定义，存储为变量Hk，Mk，Qk。
所述的数据准备与存储具体为：
21)自行根据参数将①中参数为“是”的推力设备进行编号为1，2，…，pr，将①中参数为“否”的推力设备进行编号为pr+1，…，p，同时将各推力设备变量与编号对应；
22)计算并存储推力系数矩阵B：B=[Br，Bf]
其中：
Br=10···1001···01-l1,yl1,x···-lpr,ylpr,xBf=cosαpr+1···cosαpsinαpr+1···sinαplpr+1···lp]]>
其中，

23)计算并存储4k组变量，mki，nki，其中k=1，2，…，pr，i＝1，2，3，4：
Dk={α|αkmin≤α≤αkmax}，k=1，2，…，pr，
Ft={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>
Dki=Dk∩Fi，
mki＝inf{tanα|α∈Dki}，nki=sup{tanα|α∈Dki}；
24)计算并存储组二次规划变量，每组编号为t，其中将t进行4进制转换：
t=(dprdpr-1…d1)4
将4进制的第(j-1)位存储为dj，其中j=1，2，…，pr；
第t组需计算和存储的约束变量为：

其中：
At＝[A1，A2，t，A3，t]T，Bt=[b1，0]T，Ct=[0，C3]T
A1=[E，-E]T，b1=[b11，b12]T
b11=[22T1,...,22Tpr,T(pr+1)max,...,Tpmax]T,b12=-b11;]]>

A2,t=A2,t,1...A2,t,prA3,t=A3,t,1...A3,t,prC3=C3,1...C3,pr]]>
A2temp=(mkie2k-1-e2k)T(-nkie2k-1+e2k)TA2,t,k=A2temp,dk=0,3-A2temp,dk=1,2C3,k=-e2kTe2kT]]>
A3temp1=(-ϵke2k-1-e2k)T(-ϵke2k-1+e2k)TA3temp2=(ϵke2k-1-e2k)T(ϵke2k-1+e2k)TA3,t,k=A3temp1,dk=0,3A3temp2,dk=1,2]]>
其中ej表示单位矩阵E的第j列，
25)计算并存储优化参数
若优化参数选项为“默认”，则选择默认的Hdefault，Mdefault，Qdefault数值赋值给Hk，Nk，Qk；若优化参数选项为“自定义”，则按照推力设备的输入信息记录Hk，Mk，Qk。
然后存储如下的优化参数矩阵：
H=diag[H1,H1,...,Hpr,Hpr,Hpr+1,Hpr+2,...,Hp]]]>
M=diag[M1,M1,...,Mpr,Mpr,Mpr+1,Mpr+2,...,Mp]]]>
Q=diag[Q1,Q1,...,Qpr,Qpr,Qpr+1,Qpr+2,...,Qp]]]>
。
所述的系统对准与优化建模具体为：
31)获取推力设备的状态数据，包括每个推力设备当前的输出推力TRk以及当前时刻推力输出的方向αRk；
32)计算u0并存储，在系统运行的过程中需要不断刷新：
u0=[TR0cosαR0,TR0sinαR0,...,TRprcosαRpr,TRprsinαRpr,TR(pr+1),TR(pr+2),...,TRp]]]>
33)输入需要的总推力数据τc，其中τc=[Xc，Yc，Nc]T
其中Xc，Yc分别为推力系统在进退、横移方向所需推进力，Nc为艏摇方向所需推进力矩；
34)利用给定数据生成个子可行域上的凸二次规划问题，第t个二次规划问题的标准形式为：
minJt(u，s)=uTHu+(u-u0)TM(u-u0)+sTQs
s.t.    τc＝Bu+s
Atu≤Bt+Ctu0
其中为松弛变量，表示控制输出和给定值之间的差值；
所述的子可行域上的凸二次规划的角度约束范围为0～90°，90°～180°，180°～270°，270°～360°的其中之一。
所述的优化问题求解具体为：
利用在线二次规划算法或离线二次规划算法，根据参数生成的个二次规划问题，每个问题会产生最优解ut，st，以及相应的指标值Jt；
比较这个二次规划的结果，选取：
tm=argmint∈IJt,I={1,2,...,4pr}]]>
并将作为最终的决策向量存储。
所述的推力分配决策输出具体为：
根据决策向量计算每个推进设备的角度和推力数值：
Tk=utmT(E2k-1+E2k)utm,k=1,2,...,prepr+kutm,k=pr+1,...,n]]>
其中ej表示单位矩阵E的第j列；Ej为正交投影矩阵，其中单位矩阵
推进设备角度计算，tm转换为4进制表示：
tm=(dprdpr-1…d1)4
根据条件：
tanαk=YkXk=e2kTutme2k-1Tutmαk∈Fi=dk+1]]>
其中Fi={α|(i-1)π2≤α≤π2i},i=1,2,3,4,]]>求解出αk的数值；
最终给出的动力分配方案为：
对于可旋转的推进器的动力分配决策为(Tk，αk)，k=1，2，…，pr；对于固定的推进器的动力分配决策为Tk，k=pr+1，…，n。
某动力定位对象有3个推进装置：其中2个为带舵的螺旋桨，1个为固定的槽道推进器；推力装置在船体坐标下的坐标为(-30，10)，(30，10)，(20，10)；带舵的螺旋桨推力装置产生推力的角度范围分别为[150°，210°]，[150°，210°]；槽道推进器的安装角度为90°，需要寻找尽可能减小能耗和保护设备的推力分配方案。
用户按照图4输入模型参数(有经验的用户可进一步修改默认参数后)，确定无误后，单击图4界面中的运行按钮，即进行优化计算，计算出适用于该动力定位对象的推力分配方案。