超空泡是通过安装在水下航行体头部的空化器生成的单体空泡，可以将航行体包裹起来实现与水分离^[1-3]。沾湿面积的降低可以减少航行体90%的摩擦阻力，使航行体的运行速度达到100 m/s量级^[4-6]。与完全浸湿的航行体不同，由于缺少与水的接触，超空泡航行体丧失了大部分浮力^[7-9]。因此，航行体重力需要通过空化器和舵面产生的升力来平衡，这对执行机构的性能提出了严格要求。另外，在重力或者外界扰动等作用之下，航行体与空泡壁发生碰撞产生复杂的非线性滑行力会直接影响航行体的稳定运动，这些都给超空泡航行体的稳定控制带来挑战^[10-12]。

Dzielski等提出了一个超空泡航行体模型^[13]，该模型既能充分反映被控对象主要特性又相对简洁，因此作为基准模型被广泛引用。文献^[14-21]等均是在Dzielski模型基础上进行的非线性动力学分析与控制研究。在超空泡航行体稳定控制方面，Dzielski较早提出了线性状态反馈方法，虽然能实现超空泡航行体的稳定控制，但控制效果并不理想，尾拍现象明显。Lin等^[16]在线性反馈控制基础上设计了切换控制策略，在无滑行力阶段采用线性反馈，在有滑行力阶段采用反馈线性化将滑行力进行补偿。Mao等^[17]和Zhao等^[18]研究了鲁棒性较强的超空泡航行体变结构控制，由于变结构控制本质上是一种不连续控制，容易出现抖振问题导致控制性能下降。陈超倩等^[19]和王京华等^[20]也利用平衡点处线性化方法分别设计了最优控制器和预测控制器。这种线性化方法控制效果理想但对模型依赖严重，在实际中存在一定局限性。Li等通过径向基函数神经网络对模型中不确定项进行逼近并基于Lyapunov稳定性理论自适应调整网络权重，保证了系统收敛性^[21]。

圆判据继承了经典控制论中Nyquist判据思想，同时可以处理多种非线性特性，能够从理论层面给出确切的稳定边界^[22]，因此在控制系统稳定性分析以及饱和吸引域估计等领域被广泛应用^[23-26]。范辉等^[27]和韩云涛等^[28]基于圆判据定理分析了超空泡航行体绝对稳定控制鲁棒性问题。

本文结合超空泡航行体模型特性设计出内-外环串级控制结构，推导了串级误差状态方程，基于圆判据定理分析了具有非线性滑行力的内环绝对稳定性问题。最后利用Nelder-Mead算法对内环反馈参数进一步优化，并通过仿真实验验证该控制方法的优越性。

1 超空泡航行体数学模型及圆判据定理

1.1 超空泡航行体数学模型

Dzielski模型的状态空间表达式为：

其中：x=[z θ w q]^T，C=[0 0 1 0]，u=[δ_e δ_c]^T，$\mathit{A}=\left[\begin{array}{cccc}0& -V& 1& 0\\ 0& 0& 0& 1\\ 0& 0& a_{33}& a_{34}\\ 0& 0& a_{43}& a_{44}\end{array}\right]，\mathit{B}=\left[\begin{array}{cc}0& 0\\ 0& 0\\ b_{31}& b_{32}\\ b_{41}& b_{42}\end{array}\right]$，E=[0 0 g 0]^T，P=[0 0 p₃ p₄]^T；z为深度，θ为俯仰角，w为纵向速度，q为俯仰角速度，δ_e和δ_c分别表示尾舵偏角和空化器偏角，V表示沿弹体轴方向速度，F_p表示滑行力；系数矩阵元素的表达式为 $a_{33}=\frac{C_{1}VT}{m}\left(\frac{1-n}{L}S+\frac{17}{36}nL\right)$， $a_{34}=SVT\left(\frac{-n}{m}{C}_1+\frac{7}{9}\right)-VT\left(\frac{-n}{m}{C}_1+\frac{17}{36}\right)\times \frac{17}{36}{L}^2$，$a_{43}=\frac{C_{1}VT}{m}\left(-\frac{17}{36}-\frac{11}{36}n\right)，a_{44}=\frac{-11C_{1}VTnL}{36m}，b_{31}=\frac{C_1{V}^{2}Tn}{m}\left(\frac{-S}{L}+\frac{17}{36}L\right)，b_{32}=\frac{C_1{V}^{2}TS}{mL}，b_{41}=-\frac{11C_1{V}^{2}Tn}{36m}，b_{42}=-\frac{17C_1{V}^{2}T}{36m}，p_3=\frac{T}{m}\left(\frac{S}{L}-\frac{17}{36}L\right)，p_4=\frac{11T}{36m};$ $S=\frac{11}{60}{R}^2+\frac{133}{405}{L}^2，T=\frac{1}{（7/9）S-（289/1296）L^2}，C_x=C_{x0}（1+\sigma ），C_1=\frac{C_x{R}_{\mathrm{n}}^2}{2R^2}$。其中V=75 m/s，重力加速度g=9.81 m/s²，航行体与周围水密度之比m=2，航行体长度L=1.8 m，空化数σ=0.03，航行体半径R=0.050 8 m，空化器半径R_n=0.019 1 m，舵与空化器有效比n=0.5，升力系数C_x0=0.82。滑行力F_p为纵向速度w的函数，即：

式中：α′为航行体尾部与空泡接触形成的浸入角；R_c为空泡与航行体碰撞位置的半径，简称空泡半径；h′为航行体浸入水中深度与航行体半径之比，称为相对浸入深度。相对浸入深度与浸入角是影响滑行力的关键物理量。

$\begin{array}{c}{h}^{\mathrm{\text{'}}}=\left\{\begin{array}{cc}0& \left|w\right|⩽\frac{V\left(R_{\mathrm{c}}-R\right)}{L}\\ \frac{L}{R}\left|\frac{w}{V}\right|-\frac{R_{\mathrm{c}}-R}{R}& \left|w\right|>\frac{V\left(R_{\mathrm{c}}-R\right)}{L}\end{array}\right.\\ {\alpha }^{\mathrm{\text{'}}}=\left\{\begin{array}{cc}\frac{w}{V}-\frac{{\dot{R}}_{\mathrm{c}}}{V}& w>0\\ \frac{w}{V}+\frac{{\dot{R}}_{\mathrm{c}}}{V}& w⩽0\end{array}\right.\\ \end{array}$

其中，空泡半径R_c以及空泡半径变化率${\dot{R}}_{\mathrm{c}}$的计算采用Logvinovich公式。首先定义两个常数：

$\begin{array}{c}{k}_1=\frac{L}{R_{\mathrm{n}}}{\left(\frac{1.92}{\sigma }-3\right)}^{-1}-1\\ k_2={\left[1-\left(1-\frac{4.5\sigma }{1+\sigma }\right)k_1^{(40/17)}\right]}^{1/2}\end{array}$

然后得到的空泡半径的公式：

$R_{\mathrm{c}}=R_{\mathrm{n}}{\left(0.82\frac{1+\sigma }{\sigma }\right)}^{0.5}{k}_2$

1.2 圆判据定理

为引出圆判据定理先给出以下定义^[23]：

定义1  系统结构图前向通道为线性定常系统，反馈部分为非线性静态环节，称该连接形式为Lur′e形式反馈连接，如图1所示。

定义2  无记忆函数$\psi :[0，\mathrm{\infty }）\times {\mathbb{R}}^P\to {\mathbb{R}}^P$，若有${\left[\psi \left(t，y_1\right)-K_{11}{y}_1\right]}^{\mathrm{T}}\left[K_{12}{y}_1-\psi \left(t，y_1\right)\right]$且${\tilde{K}}_1=K_{11}-K_{12}={\tilde{K}}_1^{\mathrm{T}}>0$，则称非线性函数ψ属于扇形区域[K_l1，K_l2]，即满足扇形条件。

定义3  考虑Lur′e形式反馈连接，假设外部输入r=0，系统状态方程为：

其中：${\mathit{x}}_1\in {\mathbb{R}}^n;{\mathit{u}}_1，{\mathit{y}}_1\in {\mathbb{R}}^P;\left({\mathit{A}}_1，{\mathit{B}}_1\right)$可控，（A_l，C_l）可观测；$\psi :[0，\mathrm{\infty }）\times {\mathbb{R}}^P\to {\mathbb{R}}^P$满足扇形区域条件。如果对于给定的扇形区域内的所有非线性特性，原点都是全局一致渐进稳定的，则系统是绝对稳定的。

定义4 D（α，β）是复平面中的闭圆盘，其直径是连接-（1/α）+j0和-（1/β）+j0两点的线段。

圆判据定理 ^[23]的具体内容如下：

考虑形如式（3）的系统，这里{A_l，B_l，C_l，D_l}是G（s）的一个最小实现，且ψ∈[α，β]。如果满足以下条件之一，则系统绝对稳定。

条件1：如果0＜α＜β，G（jω）的奈奎斯特曲线不进入圆盘D（α，β）内，且沿逆时针方向环绕τ次，其中τ是G（s）具有正实部的极点数。

条件2：如果0=α＜β，G（s）是赫尔维茨的，且G（jω）的奈奎斯特曲线位于直线Re[s]=-1/β右侧。

条件3：如果α＜0＜β，G（s）是赫尔维茨的，且G（jω）的奈奎斯特曲线位于圆盘D（α，β）内部。如果仅在一个区间[a，b]内满足扇形区域条件，则上述条件保证了系统在有限区域内绝对稳定。

注：考虑圆判据条件③，当β无限接近0时，圆盘D（α，β）半径会随着-（1/β）+j0而无限增大，因此当β=0时可以得到一个类似圆判据条件②的推论——如果α＜β=0，G（s）是赫尔维茨的，且G（jω）的奈奎斯特曲线位于直线Re[s]=-1/α左侧，则系统绝对稳定。

圆判据提供了一个仅通过G（jω）的奈奎斯特曲线就可以判断系统稳定性的方法，同样给出G（jω）的奈奎斯特曲线，也可确定系统绝对稳定所允许的扇形区域。

2 超空泡航行体稳定控制

根据超空泡航行体Dzielski模型，如式（1）所示，控制量可以直接作用到状态[w，q]^T，状态[w，q]^T又可以直接影响状态[z，θ]^T，因此可以将原系统转化为内-外环串级控制结构，外环通过比例控制生成虚拟指令作为内环系统的参考输入。超空泡航行体串级控制结构如图2所示，其中K₁、K₂为2×2反馈增益矩阵。

对于Dzielski模型，令η=[z，θ]^T，ξ=[w，q]^T，重新设置${\mathit{A}}_{11}=\left[\begin{array}{cc}0& -V\\ 0& 0\end{array}\right]，{\mathit{A}}_{12}=\left[\begin{array}{c}10\\ 01\end{array}\right]，{\mathit{A}}_{22}=\left[\begin{array}{c}{a}_{33}{a}_{34}\\ a_{43}{a}_{44}\end{array}\right]，{\mathit{B}}^{\mathrm{\text{'}}}=\left[\begin{array}{c}{b}_{31}{b}_{32}\\ b_{41}{b}_{42}\end{array}\right]，{\mathit{C}}^{\mathrm{\text{'}}}={\left[\begin{array}{c}1\\ 0\end{array}\right]}^{\mathrm{T}}，{\mathit{E}}^{\mathrm{\text{'}}}=\left[\begin{array}{c}g\\ 0\end{array}\right]，{\mathit{P}}^{\mathrm{\text{'}}}=\left[\begin{array}{c}{p}_3\\ p_4\end{array}\right]$。

原Dzielski模型可以转化为：

为研究跟踪问题，将式（4）转化为误差状态方程。设状态误差η_e=η-η_d，ξ_e=ξ-ξ_d，其中η_d、ξ_d分别为状态设定值。则

由式（5）可以看出，重新推导的误差状态方程并非典型的串级控制形式，为此对式（5）进行等价变换。外环采用状态反馈控制律，设ξ_d=-K₁η_e，并令

则

对χ_e求导可得：

将式（5）、式（6）和ξ_d=-K₁η_e代入式（8），并令

$\begin{array}{c}\lambda =-{\mathit{A}}_{22}{\mathit{A}}_{12}^{-1}\left({\mathit{A}}_{11}{\mathit{\eta }}_{\mathrm{d}}-{\dot{\mathit{\eta }}}_{\mathrm{d}}\right)+{\mathit{A}}_{22}{\mathit{\xi }}_{\mathrm{d}}+\\ {\mathit{K}}_1\left({\mathit{A}}_{11}-{\mathit{A}}_{12}{\mathit{K}}_1\right){\mathit{\eta }}_{\mathrm{e}}+{\mathit{A}}_{12}^{-1}\left({\mathit{A}}_{11}{\dot{\mathit{\eta }}}_{\mathrm{d}}-{\ddot{\mathit{\eta }}}_{\mathrm{d}}\right)\end{array}$

则

通过设置控制量，使u=-B′^-1（E′+λ）+υ，将常数项补偿掉，则有

式中，υ为内环误差状态方程的参考输入。为方便起见，令M=A₁₁-A₁₂K₁，N=K₁A₁₂+A₂₂，联立式（7）、式（10）可得具有典型串级控制形式的误差状态方程

一般来说，外环状态反馈矩阵K₁要保证M为赫尔维兹矩阵，选择范围比较大。另一个重要问题是如何设计内环控制律保证χ_e能收敛为0。

针对本文中Dzielski模型采用的滑行力F_p（w）分析系统的非线性反馈环节。对于输出变量y=C′χ_e，考虑到η_d=[z_d θ_d]^T=[z_d 0]^T，${\dot{\mathit{\eta }}}_{\mathrm{d}}$=[${\dot{z}}_{\mathrm{d}}$ ${\dot{\theta }}_{\mathrm{d}}$]^T=[0 0]^T，由式（11）可得y=C′χ_e=C′ξ_e。又因为C′=[1 0]，纵向速度设定值为0，因此y=w。然后，将式（11）第二个等式转化为前向通道带有状态反馈的Lur′e连接形式。设υ=-K₂χ_e，则内环Lur′e连接形式为：

代入参数值计算可得W=-B′^-1P′=[0.006 0]，所以内环系统前向通道仍为单输入单输出系统。非线性反馈环节WF_p与纵向速度的关系及其与纵向速度比值的变化曲线如图3所示。可以看出，非线性反馈环节斜率有界，因此属于扇形[α，β]，α=-0.476 2＜β=0。因此通过设置反馈矩阵K₂使前向通道传递函数是赫尔维兹的，同时满足其奈奎斯特曲线位于直线Re[s]=2.1左侧，则可保证内环系统为绝对稳定的，实现对虚拟指令的跟踪。

3 基于Nelder-Mead算法反馈参数优化

为进一步辅助内环系统参数整定使其避免陷入执行机构饱和，这里提出一种基于Nelder-Mead算法的反馈参数优化方法。Nelder-Mead算法是针对无确定解或求解复杂度高的问题所提出的一种单纯形搜索算法，该算法通过问题自身的信息进行扩展，向最有利的方向探索问题的解，具有较高搜索效率^[29]。对κ维待优化向量，选择κ+1个参数点，每个参数点是一个κ维向量（表示对待优化向量的估计）。这些参数点作为κ维向量空间中单纯形的顶点，每个点对应一个目标函数值。根据目标函数值的大小自主改变单纯形的边长和顶点，使单纯形向着目标函数值减小的方向演化，直到单纯形的各顶点对应的目标函数值无明显差异^[30]。选择二次型指标函数：

式中，Q、T均为实对称矩阵，分别表示跟踪误差和控制量输入的权重。采用Nelder-Mead搜索方法求解流程如算法1所示。

4 仿真验证与对比分析

目前研究中已有比较成熟的方法可以将超空泡航行体基准模型直接转化为Lur′e反馈形式，如文献^[27-28]等先将式（1）转化为误差状态方程，然后通过设计控制量u，实现了对系统的Lur′e反馈形式转变，进而分析超空泡航行体的绝对稳定性问题。带有状态反馈的超空泡航行体Lur′e形式如图4所示，其中K为反馈矩阵，此处W=-（B^TB）^-1B^TP。基于本文中采用的滑行力模型，设超空泡航行体初始状态为x₀=[0 0 2 0]^T，跟踪目标为x_d=[1 0 0 0]^T。设置三组仿真S1、S2、S3，选择期望的极点分别为0.5X、0.7X、X，其中X=[-30 -45 -50 -60]。传统非串级结构下分析超空泡航行体绝对稳定性的仿真结果如图5、图6所示。由图5可以看出，随着配置极点远离虚轴，带有状态反馈的前向传递函数的奈奎斯特曲线的最右边越来越接近全局绝对稳定的边界。同时，航行体状态变化曲线逐渐由发散、等幅振荡过渡到收敛，如图6所示。反馈矩阵K对原系统极点进行配置，实现了超空泡航行体的绝对稳定控制，符合圆判据定理。然而，随着系统稳定性加强，控制量输出也会逐渐增大并超出合理范围，俯仰角变化也过于剧烈，如图5（d）所示，在5 ms内便接近200（°）/s的量级，这种绝对稳定性已经失去了现实意义。

算法1 反馈参数优化算法

Alg.1 Feedback parameter optimization algorithm

超空泡航行体串级控制结构具有外环和内环两个反馈矩阵K₁与K₂，参数调节更加灵活。设超空泡航行体初始状态为x₀=[0 0 3 0]^T，跟踪目标为x_d=[1 0 0 0]^T。设置四组S4、S5、S6、S7仿真，选择反馈矩阵K₁使得外环系统矩阵M期望极点分别为1.5×[-5 -8]，1.2×[-5 -8]，0.9×[-5 -8]，0.6×[-5 -8]。内环反馈矩阵K₂=-$\left[\begin{array}{cc}-0.15& -0.3\\ 0.2& -0.2\end{array}\right]$保持固定不变。内环Lur′e形式前向传递函数奈奎斯特曲线如图7所示，由图7可以看出，外环反馈矩阵K₁对内环Lur′e系统前向传递函数有影响，但不改变内环Lur′e系统的绝对稳定性。图8为随外环反馈矩阵变化的滑行力曲线，图9为相应航行体状态变化曲线，图10为航行体执行器偏转曲线。可以看出，随着配置极点远离虚轴，跟踪调节时间越短，控制量输出越大，反之调节时间越长，控制量输入越小。因此串级结构控制下反馈矩阵K₁能够在内环绝对稳定的情况下调节过渡时间、避免执行机构饱和。本文中，设执行机构偏角允许范围为±30°，由图10可知，在仿真S5、S6、S7中未出现执行机构饱和。

基于Nelder-Mead法的内环控制参数优化方法可以进一步作为辅助调参手段避免执行机构陷入饱和。在仿真S7基础上，将跟踪目标设置为x_d=[2.5 0 0 0]^T，保持反馈矩阵K₁不变，以${\mathit{K}}_2=\left[\begin{array}{cc}-0.15& -0.3\\ 0.2& -0.2\end{array}\right]$作为初值进行内环控制参数优化，选择$\mathit{Q}=\left[\begin{array}{cc}20& 0\\ 0& 20\end{array}\right]，\mathit{T}=\left[\begin{array}{cc}200& 0\\ 0& 200\end{array}\right]$，单纯形边长a=0.1，结束搜索条件ε=0.01。经过24次迭代，最终确定内环控制参数${\mathit{K}}_2=\left[\begin{array}{cc}-0.0426& -0.2370\\ 0.0768& -0.4029\end{array}\right]$，仿真结果如图11~图13所示。可以看出，优化前后内环系统都是绝对稳定的，航行体最终实现深度的稳定跟踪，但优化前的航行体在运动过程中出现控制量饱和。由图13可知，优化前空化器偏转最大40°左右，尾舵偏转最大-50°左右。经过优化之后确定的控制参数使得控制输出明显减小，均处于执行机构允许范围内。Nelder-Mead算法目标函数值变化曲线如图14所示，由图14可以看出，随着每次迭代的进行，目标函数值向着减小的方向变化。

5 结论

本文基于圆判据理论和Nelder-Mead算法提出了超空泡航行体深度跟踪串级控制方法。根据Dzielski模型特性将原系统转化为内-外环串级控制结构，使得参数调节更加灵活。内环系统可以根据圆判据定理保证其绝对稳定性，并利用Nelder-Mead算法对内环反馈参数进一步优化，充分利用有限控制量实现对深度的稳定跟踪。对比仿真结果验证了本文所提方法的有效性。由于采用串级控制结构，为保证整个系统稳定，内环系统要求具有较高带宽。在将来的研究工作中，将重点研究超空泡航行体串级控制结构下的全局稳定性问题。

参考文献