Optimal Scheduling for Regional Integrated Energy System Operation Based on the AMOWOA
-
摘要: 目前, 智能优化算法已广泛应用于工程优化中, 在当前多能耦合与互补的能源发展趋势下, 仅考虑系统经济指标的单目标优化模式已经不再适用于目前区域综合能源系统(Integrated energy system, IES)的运行优化调度, 需要研究一种多目标运行策略来解决区域综合能源系统的运行优化调度问题. 首先综合考虑经济与能源利用两个指标并结合商业住宅区域的特性, 以系统日运行收益和一次能源利用率为优化目标构建商业住宅区域综合能源系统多目标运行优化调度模型. 其次由于传统多目标智能优化算法缺乏一种最优解综合评价方法, 基于非支配排序以及拥挤度计算的多目标算法框架, 提出一种利用模糊一致矩阵选取全局最优解的多目标鲸鱼优化算法(A multi-objective whale optimization algorithm, AMOWOA), 并将提出算法对商住区域综合能源系统多目标运行优化调度模型进行求解. 最后以华东某商业住宅区域综合能源系统为例进行仿真, 验证了该方法的有效性和可行性.Abstract: Now intelligent optimization algorithms have been widely used in engineering optimization. Under the current energy development trend of multi-energy coupling and complementation, the single-objective optimization model that only considers system economic indicators is no longer applicable to the optimal scheduling for integrated energy system (IES) operation, it is necessary to study a multi-objective operation strategy to solve operation optimization and scheduling of regional integrated energy system problem. First, considering the two indicators of economy and energy utilization and combining the characteristics of commercial and residential areas, a multi-objective operation optimization and scheduling model of an integrated energy system for commercial and residential areas with system daily operating income and primary energy utilization as the optimization goals is constructed. Secondly, since the traditional multi-objective intelligent optimization algorithm lacks a comprehensive evaluation method of optimal solutions, based on the framework of multi-objective algorithm with non-dominated sorting and congestion calculation, a multi-objective whale optimization algorithm (AMOWOA) is proposed to select the optimal solution using fuzzy consistency matrix, and then the optimal solution will be evaluated by AMOWOA, and then the optimal solution will be selected by AMOWOA. AMOWOA is proposed to select the global optimal solution using a fuzzy consistency matrix, and the proposed algorithm is used to solve the multi-objective operation optimisation scheduling model of an integrated energysystem in a commercial and residential area. The proposed algorithm is used to solve the multi-objective operation optimisation model of commercial and residential energy system. Finally, a commercial and residential area integrated energy system in East China is used as an example to verify the effectiveness and feasibility of the proposed method.
-
高超声速变外形飞行器(Hypersonic morphing vehicle, HMV)是一种环境适应能力和突防生存能力较强的新型可变构型高超声速飞行器. 飞行器飞行过程中, 变形会引起气动力、气动力矩和质心、压心、转动惯量等参数的改变, 导致飞行器存在较大的模型不确定性; 同时, 飞行过程中外界环境复杂且变化剧烈, 飞行器不仅受到外界未知干扰的影响, 还需要满足面向特定任务与复杂环境的性能约束条件. 因此, 设计高超声速变外形飞行器高精度、强鲁棒和满足指定性能约束要求的姿态控制方法具有重要意义.
通常可将飞行器模型不确定性和飞行过程中的外部干扰视为复合总扰动, 再通过扰动观测的方法对其进行精确估计与补偿, 从而提高控制系统的鲁棒性. 基于扩张状态观测器[1-2]、干扰观测器[3-4]的控制方法在现有变外形飞行器姿态控制研究中获得了广泛应用. 文献[1]将折叠翼飞行器纵向非线性动力学模型中存在的非线性项、耦合项以及参数时变项都视为系统内外总扰动, 利用扩张状态观测器对总扰动进行实时估计和补偿, 实现了飞行器的高精度稳定控制. 文献[3]基于模糊理论提出了一种模糊干扰观测器, 可在有限时间内对飞行器模型的匹配和非匹配扰动进行精确估计, 求解得到前馈控制补偿用于变外形飞行器鲁棒姿态控制. 文献[4]建立了固定时间干扰观测器以实现对变外形飞行器包含模型不确定性、外部干扰和执行机构故障信息的复合干扰的快速精确估计, 并基于反步法设计了具有较强抗干扰能力的容错控制器.
飞行器控制系统在满足稳定性的前提下往往还存在快速性和高精度的性能要求[5]. 预设性能控制方法因其在定量表征控制系统的瞬态性能和稳态性能方面的突出优势, 被广泛应用于高超声速飞行器控制[6-7]. 文献[8]针对干扰有界但未知的可重复使用飞行器姿态控制系统, 提出了一种基于预设性能的自适应多变量控制算法, 可以保证再入姿态跟踪误差被限制在与有界扰动无关的预设区域. 文献[9]综合考虑可重复使用飞行器执行器饱和、气动参数摄动和外部扰动, 融合高增益扩张状态观测器和低复杂度的输出反馈扰动补偿控制方案, 提出了一种基于反步法的保性能姿态控制方法. 文献[10]研究了存在跟踪误差约束和模型不确定性影响下的高超声速飞行器控制问题, 提出一种新型预设时间性能函数以确保跟踪误差在预定义时间内收敛至预设约束边界, 实现了飞行器快速且高精度的跟踪控制.
目前针对变外形飞行器变形及控制的研究多集中于飞行速度小、高度低的航空飞行器[11-12]. 文献[13]研究了非匹配扰动影响下的半倾转旋翼变外形无人机的轨迹跟踪控制问题. 文献[14]设计开发了一种新的仿生变外形无人机并实现了变外形无人机构型优化的自主变形控制. 针对可变构型高超声速飞行器, 文献[2]研究了一类可变后掠飞行器的建模与姿态控制问题, 文献[15]研究了一类可变展长飞行器的制导控制一体化问题.
总体来看, 当前针对可变构型高超声速飞行器建模及控制问题的研究还相对较少, 尤其是采用大尺度折叠变构方式的高超声速变外形飞行器, 鲜有公开资料. 因此, 开展此类飞行器变形影响下具有强耦合、强非线性特征的运动建模和气动特性分析具有重要现实意义. 同时, 飞行器飞行过程中存在的模型不确定性和外界干扰具有未知、复杂、多源以及与飞行状态和变外形过程强烈耦合等特点, 这对所设计控制方案的干扰适应能力提出了较高要求. 因此, 可采用收敛时间具有常值上界的有限时间观测器[16]或固定时间观测器[17]对扰动进行精确估计和前馈补偿, 既为提高控制系统鲁棒性提供了简单、高效的解决方案, 也能降低观测器调节过程对系统的不利影响. 进一步, 考虑到飞行器飞行过程中的特定任务和复杂环境, 预设性能控制为确保姿态控制系统具备期望的瞬态及稳态性能提供了良好的解决途径[18]. 因此, 可基于预设性能控制框架设计姿态控制器, 实现姿态跟踪误差在预定义时间内收敛至平衡点的给定邻域, 并确保在不同变形条件及干扰影响下均能实现期望的控制性能.
本文针对大尺度变外形、模型不确定性和外部干扰影响下的高超声速变外形飞行器建模及控制问题开展研究, 提出了一种可实现预先设定瞬态和稳态性能的高精度、强鲁棒固定时间控制方法. 本文的主要贡献如下:
1) 本文研究对象为高超声速下的变外形滑翔飞行器, 具有一定的特殊性, 且目前对该研究对象的具体建模与控制的研究还较少, 缺乏相应的参考资料. 鉴于此, 本文在该飞行器的建模及特性分析方面进行了较为完备的工作, 具体包括: 建立了变形过程中由机翼折叠引起的附加力及附加力矩的具体形式和折叠翼变质心的位置方程; 建立了飞行器气动模型并开展气动特性分析, 分析了变折叠角对飞行器气动特性的影响规律, 研究了机翼折叠变形以及附加力/附加力矩对飞行器闭环动态特性的影响特性. 本部分工作不仅是本文姿态控制方法设计的重要基础, 也可为高超声速变外形滑翔飞行器轨迹规划、制导与控制提供较为可靠的模型依据.
2) 采用固定时间干扰观测器来估计模型不确定性和外部干扰构成的复合总扰动, 提升了系统的抗干扰能力. 提出了一种新型固定时间预设性能函数, 定量描述姿态跟踪误差的瞬态和稳态性能, 并基于预设性能控制和动态面控制设计了姿态跟踪控制器, 实现系统的固定时间收敛. 与现有基于滑模控制的固定时间控制方法相比, 该方法避免了状态量的分数幂次运算和符号函数, 既减少了计算量, 也克服了由符号函数导致的控制量不连续的问题. 同时, 该方法采用了积分障碍Lyapunov函数处理性能约束, 并将其嵌入动态面控制和预设性能控制相结合的控制方案中, 既简化了设计流程, 又保证了期望的优异控制性能.
本文其余内容安排如下: 第1节建立了高超声速变外形飞行器的运动模型、气动模型和姿态控制模型; 第2节提出了一种基于预设性能的固定时间姿态控制方法, 包括固定时间干扰观测器、预设性能函数和基于动态面控制的预设性能控制器的设计; 第3节进行了仿真验证和结果分析; 第4节总结了本文工作并得出结论.
1. 飞行器运动建模及问题描述
1.1 研究对象概述
本文以图1所示的折叠式高超声速变外形飞行器为研究对象. 作为一种无动力飞行器, 其通常由助推器运送到某一高度并分离, 依靠自身较大的升阻比, 借助气动升力以实现宽速域、大空域和长航程的高速滑翔飞行.
图 1 高超声速变外形飞行器气动外形及变形过程示意图Fig. 1 Aerodynamic shape and morphing process of hypersonic morphing vehicle作为一种含有大尺度变形机构的飞行器, 高超声速变外形飞行器可根据任务需要和环境变化自主改变气动外形, 以保证飞行器在任意飞行阶段均具有优良的气动特性和飞行性能. 折叠式高超声速变外形飞行器是其中的一种典型形式, 能够通过机翼的折叠改变自身构型, 以适应于不同的飞行任务. 机翼折叠时, 翼面积减小, 可有效减小阻力, 适合高速飞行及大范围机动; 机翼展开时, 翼面积增大, 升阻比增大, 适合高空长航时滑翔, 可有效提高射程.
相较于传统固定构型飞行器, 折叠式高超声速变外形飞行器在机身前部多出一对可变形的折叠翼, 变形过程如图1所示. 其中, 折叠翼1和折叠翼2同时运动且对称变形, 当$ {{\delta}_{f}} = {{0}^{\circ}} $时, 机翼为完全展开状态; 当$ {{\delta}_{f}} = {{155}^{\circ}} $时, 机翼为完全收起状态.
本文仅针对30 ~ 40 km高度的无动力滑翔段开展研究, 重点研究该阶段下考虑机翼折叠变形、模型不确定性和外部干扰影响的高性能姿态控制方法.
1.2 飞行器运动模型
考虑如下假设:
假设 1[19]. 忽略地球自转, 将地球视为均质圆球, 大气相对地球静止且在相同高度均匀分布.
假设 2[19]. 忽略姿态运动方程中的质心运动耦合项, 即, 认为$ \dot{\lambda} = \dot{\phi} = \dot{\theta} = {{\dot{\psi}}_{v}} = 0 $.
假设 3[1-2]. 将飞行器机体视作主刚体, 两侧折叠翼视作从刚体, 整个飞行器系统可视为由3个刚体组成的多刚体系统. 飞行器折叠翼做对称变形, 且翼密度及厚度均匀.
根据假设1、2, 高超声速变外形飞行器六自由度运动方程可表示如下:
$$ \left\{\begin{aligned} &\dot{R} = {} V\sin\theta \\ & \dot{\lambda} = {} -\frac{V\cos\theta\sin{{\psi}_{v}}}{R\cos\phi} \\ & \dot{\phi} = {} \frac{V\cos\theta\cos{{\psi}_{v}}}{R} \\ &\dot{V} = {} \frac{{{X}_{t}}}{m}-g\sin\theta \\ & \dot{\theta} = {} \frac{{{Y}_{t}}\cos\sigma-{{Z}_{t}}\sin\sigma}{mV}- \frac{g\cos\theta}{V}+\frac{V\cos\theta}{R} \\ & {{{\dot{\psi}}}_{v}} = {} -\frac{{{Y}_{t}}\sin\sigma+{{Z}_{t}}\cos\sigma}{mV\cos\theta}+ \frac{V\cos\theta\sin{{\psi}_{v}}\tan\phi}{R} \end{aligned}\right. $$ (1) $$ \left\{\begin{aligned} &\dot{\alpha} = {} \omega_z-\tan\beta\left(\omega_x\cos\alpha-\omega_y\sin\alpha\right) \\ &\dot{\beta} = {} \omega_x\sin\alpha+\omega_y\cos\alpha \\ &\dot{\sigma} = {} \sec\beta\left(\omega_x\cos\alpha-\omega_y\sin\alpha\right) \\ & \dot{\omega}_x = {} I_2 M_{tx}+I_4 M_{ty}+I_6\omega_x\omega_z+I_7\omega_y\omega_z\\ &\dot{\omega}_y = {} I_4 M_{tx}+I_1 M_{ty}+I_8\omega_x\omega_z-I_6\omega_y\omega_z\\ & \dot{\omega}_z = {} I_3 M_{tz}+I_5 \omega_x^2-I_5\omega_y^2+I_9\omega_x\omega_y \end{aligned}\right. $$ (2) 式中, $ R $为地心距; $ m $, $ V $分别为飞行器质量和速度; $ \lambda $, $ \phi $分别为经度和纬度; $ \theta $, $ {{\psi }_{v}} $分别为航迹倾角和航迹偏角; $ \alpha $, $ \beta $, $ \sigma $分别为攻角、侧滑角和倾侧角; $ {{\omega }_{x}} $, $ {{\omega }_{y}} $, $ {{\omega }_{z}} $为机体轴三通道角速率; $ {{I}_{1}} $ ~ $ {{I}_{9}} $为飞行器惯量参数, 具体形式见式 (3); $ {{X}_{t}} $, $ {{Y}_{t}} $, $ {{Z}_{t}} $分别为总阻力、总升力和总侧力, $ {{M}_{tx}} $, $ {{M}_{ty}} $, $ {{M}_{tz}} $分别为滚转、偏航和俯仰通道的合力矩, 其具体形式见式 (4).
$$ \left\{\begin{aligned} I_*& = I_{xx}I_{yy}-I_{xy}^2 , I_1 = \frac{I_{xx}}{I_*},\\ I_2& = \frac{I_{yy}}{I_*}, I_3 = \frac{1}{I_{zz}} , I_4 = \frac{I_{xy}}{I_*}, I_5 =\frac{ I_{xy}}{I_{zz}} \\ I_6& = I_{xy}\frac{I_{zz}-I_{xx}-I_{yy}}{I_*} \\ I_7& = \frac{I_{xy}^2+I_{yy}^2-I_{yy}I_{zz}}{I_* }\\ I_8& = \frac{I_{xx}I_{zz}-I_{xx}^2-I_{xy}^2}{I_*} , I_9 = \frac{I_{xx}-I_{yy}}{I_{zz}} \end{aligned}\right. $$ (3) $$ \left\{ \begin{aligned} & \underbrace{\left[ \begin{matrix} {{X}_{t}} \\ {{Y}_{t}} \\ {{Z}_{t}} \\ \end{matrix} \right]}_{{{\boldsymbol{F}}_{t}}}=\underbrace{\left[ \begin{matrix} -D \\ L \\ Y \\ \end{matrix} \right]}_{{{\boldsymbol{F}}_{a}}}+\underbrace{\left[ \begin{matrix} {{F}_{sx}} \\ {{F}_{sy}} \\ {{F}_{sz}} \\ \end{matrix} \right]}_{{{\boldsymbol{F}}_{s}}} =\left[ \begin{matrix} {-{C_D}{Q_A}{S_r}+F_{sx}} \\ {{C_L}{Q_A}{S_r}+F_{sy}} \\ {{C_Y}{Q_A}{S_r}+F_{sz}} \\ \end{matrix} \right]\\ \\ & \underbrace{\left[ \begin{matrix} {{M}_{tx}} \\ {{M}_{ty}} \\ {{M}_{tz}} \\ \end{matrix} \right]}_{{{\boldsymbol{M}}_{t}}} = \underbrace{\left[ \begin{matrix} {{M}_{x}} \\ {{M}_{y}} \\ {{M}_{z}} \\ \end{matrix} \right]}_{{{\boldsymbol{M}}_{a}}}+\underbrace{\left[ \begin{matrix} {{M}_{sx}} \\ {{M}_{sy}} \\ {{M}_{sz}} \\ \end{matrix} \right]}_{{{\boldsymbol{M}}_{s}}} = \left[ \begin{matrix} {{C_{mx}}{Q_A}{S_r}{b_A}+M_{sx}} \\ {{C_{my}}{Q_A}{S_r}{b_A}+M_{sy}} \\ {{C_{mz}}{Q_A}{S_r}{c_A}+M_{sz}} \\ \end{matrix} \right]\\ \end{aligned} \right. $$ (4) 式中, $ D $, $ L $, $ Y $分别为阻力、升力和侧力; $ C_D $, $ C_L $, $ C_Y $分别为阻力系数、升力系数和侧力系数; $ M_x $, $ M_y $, $ M_z $分别为滚转力矩、偏航力矩和俯仰力矩; $ C_{mx} $, $ C_{my} $, $ C_{mz} $分别为滚转力矩系数、偏航力矩系数和俯仰力矩系数; $ S_r $为参考面积, $ b_A $为参考气动展长, $ c_A $为参考气动弦长; $ Q_A = {{{\rho}_{A}}{{V}^{2}}}/{2} $为动压, $ {{\rho }_{A}} $为大气密度, 所采用大气模型见文献[20]. $ {{\boldsymbol{F}}_{s}} $, $ {{\boldsymbol{M}}_{s}} $的计算公式见式(5).
$$ \left\{ \begin{aligned} & {{\boldsymbol{F}}_{s}} = {} {{\boldsymbol{T}}_{\text{VB}}}\sum\limits_{i = 1}^{2}{{{\boldsymbol{F}}_{Bsi}}},\text{ }{{\boldsymbol{M}}_{s}} = \sum\limits_{i = 1}^{2}{{{\boldsymbol{M}}_{si}}}\text{ } \\ & {{\boldsymbol{F}}_{Bsi}} = {} {{m}_{i}}\left( \frac{{{\delta }^{2}}{{\boldsymbol{s}}_{i}}}{\delta {{t}^{2}}}+2{\boldsymbol{\omega }}\times \frac{\delta {{\boldsymbol{s}}_{i}}}{\delta t}\;+ \right. \\ & \qquad\quad{} \left. \frac{\delta {\boldsymbol{\omega }}}{\delta t}\times {{\boldsymbol{s}}_{i}}+{\boldsymbol{\omega }}\times \left( {\boldsymbol{\omega }}\times {{\boldsymbol{s}}_{i}} \right) \right) \\ & {{\boldsymbol{M}}_{si}} = {}{{\boldsymbol{s}}_{i}}\times \left( {{m}_{i}}{{\boldsymbol{T}}_{\text{BG}}}{\boldsymbol{g}} \right)- {{m}_{i}}{{\boldsymbol{s}}_{i}}\times \Bigg( {{\boldsymbol{T}}_{\text{BH}}}\frac{\delta {{\boldsymbol{v}}_{0}}}{\delta t}\;+ \\ &\qquad\quad {} {\boldsymbol{\omega }}\times\left( {{\boldsymbol{T}}_{\text{BH}}}{{\boldsymbol{v}}_{0}} \right) \Bigg)-{{\boldsymbol{s}}_{i}}\times {{\boldsymbol{F}}_{Bsi}} \end{aligned} \right. $$ (5) 式中, $ i $ = 1, 2, 分别表示折叠翼$ i $; ${\boldsymbol{\omega}} = [ {{\omega }_{x}},{{\omega }_{y}}, {{\omega }_{z}} ]^{\text{T}}$为角速度向量; $ {{\boldsymbol{v}}_{0}} = {{[ V,0,0 ]}^{\text{T}}} $为速度向量; $ {{\boldsymbol{s}}_{i}} = {{[ {{s}_{xi}},{{s}_{yi}},{{s}_{zi}} ]}^{\text{T}}} $为折叠翼$ i $质心位置向量, 其具体推导过程见下文; $ {\boldsymbol{g}} = {{[ 0,-g,0 ]}^{\text{T}}} $为引力向量; $ {{\boldsymbol{T}}_{\text{VB}}} $, $ {{\boldsymbol{T}}_{\text{BG}}} $, $ {{\boldsymbol{T}}_{\text{BH}}} $均为转换矩阵, 具体形式见式 (6).
$$ \left\{\begin{aligned} {\boldsymbol{T}}_{\text{BV}}& = {\boldsymbol{M}}_{3}\left(\alpha\right){\boldsymbol{M}}_{2}\left(\beta\right) \\ {\boldsymbol{T}}_{\text{VH}}& = {\boldsymbol{M}}_{1}\left(\sigma\right) \\ {\boldsymbol{T}}_{\text{HG}}& = {\boldsymbol{M}}_{3}\left(\theta\right){\boldsymbol{M}}_{2}\left(\psi_{v}\right) \\ {\boldsymbol{T}}_{\text{VB}}& = {\boldsymbol{T}}_{\text{BV}}^{\text{T}} \\ {\boldsymbol{T}}_{\text{BH}}& = {\boldsymbol{T}}_{\text{BV}}{\boldsymbol{T}}_{\text{VH}} \\ {\boldsymbol{T}}_{\text{BG}}& = {\boldsymbol{T}}_{\text{BV}}{\boldsymbol{T}}_{\text{VH}}{\boldsymbol{T}}_{\text{HG}} \end{aligned}\right. $$ (6) 式中, $ {\boldsymbol{M}}_{1} $, $ {\boldsymbol{M}}_{2} $, $ {\boldsymbol{M}}_{3} $为初等转换矩阵, 其形式见文献[21].
飞行器两侧机翼对称, 则$ s_{x1} = s_{x2} $, $ s_{y1} = s_{y2} $, $ s_{z1} = -s_{z2} $. 以$ O_1 $为坐标系原点建立折叠翼参考坐标系$ {{\boldsymbol{E}}_1} $$ (O_1x_1y_1z_1) $, 如图2所示. $ \triangle {ACD}_1 $表示折叠翼2完全展开, $ \delta_{f_0} = 0^{\circ } $; $ \triangle {ACD}_2 $表示折叠角为$ \delta_{f} $; $ \triangle {ACB} $表示折叠翼完全收起, $ \delta_{f_{\max }} = 155^{\circ } $. 根据假设3, 当折叠角逐渐增大时($ \delta_{f_0}\to \delta_{f} \to \delta_{f_{\max }} $), 其质心位置变化为($ G_1\to G_2 \to G_n $), 且$ EG_1 = {ED_1}/3 $, $ E $为$ AC $中点. 过$ G_1 $作垂线$ G_1F\bot CA $, 过$ G_2 $作垂线$ {G_2}{G_p}\bot G_1F $, 则$ G_p $为折叠翼质心位置在$ \triangle {ACD}_1 $所构成平面上的投影点, 当折叠角变化时, 即$ \delta_f \in \left[ 0,155^{\circ } \right] $, $ G_p $由$ G_1 $点向$ F $点方向移动.
根据图2可知, $ O_1A $, $ O_1B $, $ O_1C $均为已知常量, 设$ O_1A = a_{h1} $, $ O_1B = b_{h1} $, $ O_1C = c_{h1} $, 令$ \angle CAD_1 = {\alpha }_{h1} $, $ \angle O_1CA = {\alpha }_{h2} $, 则有如下关系:
$$ \left\{ \begin{aligned} {{x}_{p}} = {}& \frac{{{a}_{h1}}}{2}-\frac{{{a}_{h1}}\cos 2{{\alpha }_{h1}}}{6}\;+\\ {}& \frac{{{a}_{h1}}\sin 2{{\alpha }_{h1}}\cos {{\alpha }_{h2}}}{6\sin {{\alpha }_{h2}}}\cos {{\delta }_{f}} \\ {{y}_{p}} = {}& \frac{{{a}_{h1}}\sin 2{{\alpha }_{h1}}}{6\sin {{\alpha }_{h2}}}\sin {{\delta }_{f}} \\ {{z}_{p}} = {}& {{c}_{h1}}-\left( \frac{{{a}_{h1}}}{2\sin {{\alpha }_{h2}}}-\frac{{{a}_{h1}}\cos 2{{\alpha }_{h1}}}{6\sin {{\alpha }_{h2}}} \right)\cos {{\alpha }_{h2}}\;+\\ {}& \frac{{{a}_{h1}}\sin 2{{\alpha }_{h1}}}{6}\cos {{\delta }_{f}} \end{aligned} \right. $$ (7) 式中, 除折叠角$ \delta_f $外, 其他参数均为常量, 故式(7) 可进一步改写为
$$ \left\{ \begin{aligned} {{s}_{x2}}& = {{p}_{x1}}+{{p}_{x2}}\cos{{\delta}_{f}} \\ {{s}_{y2}}& = {{p}_{y1}}+{{p}_{y2}}\sin{{\delta}_{f}} \\ {{s}_{z2}}& = {{p}_{z1}}+{{p}_{z2}}\cos{{\delta}_{f}} \end{aligned} \right. $$ (8) 式中, $ {{p}_{x1}} $, $ {{p}_{x2}} $, $ {{p}_{y1}} $, $ {{p}_{y2}} $, $ {{p}_{z1}} $, $ {{p}_{z2}} $均为常量, 与式(7) 中参数一一对应.
1.3 飞行器气动模型
将影响飞行器气动力系数和气动力矩系数的主要变量定义为$ {{\boldsymbol{X}}_{A}} $, 且$ {{\boldsymbol{X}}_{A}} = {{\left[ {{x}_{1}}, \cdots , {{x}_{k}}, \cdots , {{x}_{7}} \right]}^{\text{T}}} $ = $ {{\left[ \text{Ma}, \alpha , \beta , {{\delta }_{x}}, {{\delta }_{y}}, {{\delta }_{z}}, {{\delta }_{f}} \right]}^{\text{T}}} $. 在此基础上建立飞行器气动参数模型, 即建立${{\boldsymbol{C}}_{{F}}}$、${{\boldsymbol{C}}_{{M}}}$与$ {{\boldsymbol{X}}_{A}} $的函数关系, 即
$$ \left\{ \begin{aligned} &{{\boldsymbol{C}}_{F}} = {{f}_{F}}\left( {{\boldsymbol{X}}_{A}} \right) \\ &{{\boldsymbol{C}}_{M}} = {{f}_{M}}\left( {{\boldsymbol{X}}_{A}} \right) \end{aligned} \right. $$ (9) 式中, $ {{\boldsymbol{C}}_{F}} = {{\left[ {{C}_{D}},{{C}_{L}},{{C}_{Y}} \right]}^{\text{T}}} $ 为气动力系数向量; ${{\boldsymbol{C}}_{M}} = {{\left[ {{C}_{mx}},{{C}_{my}},{{C}_{mz}} \right]}^{\text{T}}}$为气动力矩系数向量; $ {{f}_{F}} $, $ {{f}_{M}} $通常为复杂非线性函数, 其数学模型一般包括代数模型、微分方程和积分模型等, 其中代数模型形式较为简单, 本文拟采用代数模型中的多项式模型, 其一般表达式为
$$ \left\{ \begin{aligned} &{{C}_{i}} = \sum\limits_{{{j}_{1}} = 0}^{{{N}_{1}}}{\cdots }\sum\limits_{{{j}_{p}} = 0}^{{{N}_{p}}}{{{c}_{i}}}(*)\prod\limits_{k = 1}^{p}{x_{k}^{{{j}_{k}}}},&\text{ }i = D,L,Y \\ &{{C}_{mi}} = \sum\limits_{{{j}_{1}} = 0}^{{{N}_{1}}}{\cdots }\sum\limits_{{{j}_{p}} = 0}^{{{N}_{p}}}{{{c}_{mi}}(*)}\prod\limits_{k = 1}^{p}{x_{k}^{{{j}_{k}}}},&\text{ }i = x,y,z\;\; \end{aligned} \right. $$ (10) 式中, $ {{j}_{k}} $为$ {{x}_{k}} $的幂次方, $ {{c}_{i}}(*) $和$ {{c}_{mi}}(*) $为气动参数对应于乘积项$\prod\nolimits_{k = 1}^{p}{x_{k}^{{{j}_{k}}}}$的气动导数; $0\le {{j}_{k}},{{N}_{p}}\le 3$, $ 1\le k,p\le 7 $. 式(10) 的详细展开式见式(11)、(12).
$$ \left\{ \begin{aligned} {{C}_{L}} = {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,L}^{0}{\boldsymbol{\alpha }}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,L}^{{{\delta }_{f1}}}{\boldsymbol{\alpha }}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,L}^{{{\delta }_{f2}}}{\boldsymbol{\alpha }} \right]{{{\boldsymbol{\delta }}}_{f}}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,L}^{{{\delta }_{x}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{x,L}}+{\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,L}^{{{\delta }_{z}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{z,L}} \\ {{C}_{D}} = {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{0}{{{\boldsymbol{\alpha }}}_{D}}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{f1}}}{{{\boldsymbol{\alpha }}}_{D}}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{f2}}}{{{\boldsymbol{\alpha }}}_{D}} \right]{{{\boldsymbol{\delta }}}_{f}}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{x}}}{{{\boldsymbol{\alpha }}}_{D}}{{{\boldsymbol{\delta }}}_{x,D}} + {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{y}}}{{{\boldsymbol{\alpha }}}_{D}}{{{\boldsymbol{\delta }}}_{y,D}}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{z1}}}{{{\boldsymbol{\alpha }}}_{D}}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,D}^{{{\delta }_{z2}}}{{{\boldsymbol{\alpha }}}_{D}} \right]{{{\boldsymbol{\delta }}}_{z,D}} \\ {{C}_{Y}} = {}&{\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,Y}^{0}{\boldsymbol{\alpha}} {\boldsymbol{\beta }}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,Y}^{{{\delta }_{f1}}}{\boldsymbol{\alpha }}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,Y}^{{{\delta }_{f2}}}{\boldsymbol{\alpha }} \right]{{{\boldsymbol{\delta }}}_{f}}{\boldsymbol{ \beta }}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,Y}^{{{\delta }_{x}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{x,Y}}+{\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,Y}^{{{\delta }_{y}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{y,Y}} \end{aligned} \right. $$ (11) $$ \left\{ \begin{aligned} {{C}_{mx}} = {}&{\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mx}^{0}{\boldsymbol{\alpha}} {\boldsymbol{\beta }}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mx}^{{{\delta }_{f1}}}{\boldsymbol{\alpha }}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mx}^{{{\delta }_{f2}}}{\boldsymbol{\alpha }} \right]{{{\boldsymbol{\delta }}}_{f}}{\boldsymbol{ \beta }}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mx}^{{{\delta }_{x}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{x,mx}} \\ {{C}_{my}} = {}&{\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,my}^{0}{\boldsymbol{\alpha}} {\boldsymbol{\beta }}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,my}^{{{\delta }_{f1}}}{\boldsymbol{\alpha }}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,my}^{{{\delta }_{f2}}}{\boldsymbol{\alpha }} \right]{{{\boldsymbol{\delta }}}_{f}}{\boldsymbol{ \beta }}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,my}^{{{\delta }_{y}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{y,my}} \\ {{C}_{mz}} = {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mz}^{0}{\boldsymbol{\alpha }}\;+ \\ {}& \left[ {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mz}^{{{\delta }_{f1}}}{\boldsymbol{\alpha }}, {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mz}^{{{\delta }_{f2}}}{\boldsymbol{\alpha }} \right]{{{\boldsymbol{\delta }}}_{f}}\;+ \\ {}& {\boldsymbol{M}}{{\boldsymbol{a}}^{\text{T}}}{\boldsymbol{C}}_{k,mz}^{{{\delta }_{z}}}{\boldsymbol{\alpha }}{{{\boldsymbol{\delta }}}_{z,mz}} \end{aligned} \right. $$ (12) 式中,
$$ \left\{ \begin{aligned} & {\boldsymbol{Ma}} = {{\left[ 1,{Ma} \right]}^{\text{T}}},\text{ }{{{\boldsymbol{\delta }}}_{f}} = {{\left[ 1,{{\delta }_{f}} \right]}^{\text{T}}}\text{, }{\boldsymbol{\alpha }} = {{\left[ 1,\alpha \right]}^{\text{T}}} \\ &{{{\boldsymbol{\alpha }}}_{D}} = {{\left[ 1,{{\alpha }^{2}} \right]}^{\text{T}}},{{\boldsymbol{\beta}} } = \left[ \beta \right] \\ & {{{\boldsymbol{\delta }}}_{x,L}} = \left[ \delta _{x}^{2} \right],\text{ }{{{\boldsymbol{\delta }}}_{z,L}} = \left[ {{\delta }_{z}} \right],\text{ }{{{\boldsymbol{\delta }}}_{x,D}} = \left[ \delta _{x}^{2} \right] \\ &{{{\boldsymbol{\delta }}}_{y,D}} = \left[ \delta _{y}^{2} \right],\text{ }{{{\boldsymbol{\delta }}}_{z,D}} = {{\left[ {{\delta }_{z}},\delta _{z}^{2} \right]}^{\text{T}}} \\ & {{{\boldsymbol{\delta }}}_{x,Y}} = \left[ {{\delta }_{x}} \right],\text{ }{{{\boldsymbol{\delta }}}_{y,Y}} = \left[ {{\delta }_{y}} \right],\text{ }{{{\boldsymbol{\delta }}}_{x,mx}} = \left[ {{\delta }_{x}} \right] \\ &{{{\boldsymbol{\delta }}}_{y,my}} = \left[ {{\delta }_{y}} \right],\text{ }{{{\boldsymbol{\delta }}}_{z,mz}} = \left[ {{\delta }_{z}} \right] \end{aligned} \right. $$ (13) 高超声速变外形飞行器的气动特性有别于传统固定构型飞行器, 对其进行气动特性分析有助于深刻理解飞行器运动特性, 为高精度姿态控制提供可靠模型基础. 基于式(10)所描述的飞行器气动模型, 可计算得到不同工况及不同变形条件下的气动力系数和气动力矩系数. 式(10)气动系数所涉及状态量范围如表1所示. 对于该飞行器而言, 机翼折叠变形改变翼面积, 继而影响其升阻特性, 因此, 有必要研究飞行器升阻特性变化规律. 面向本文所研究的姿态控制问题, 还应分析飞行器气动力矩系数的变化规律.
表 1 气动模型状态量范围Table 1 State quantity range of aerodynamics model状态量 符号 取值范围 马赫数 Ma $\left [ 2,18 \right ]$ 攻角 $\alpha $ $\left [ 0^{\circ},20^{\circ} \right ]$ 侧滑角 $\beta $ $\left [ -2^{\circ},2^{\circ} \right ] $ 滚转舵偏角 $\delta_x$ $ \left [ -20^{\circ},20^{\circ} \right ] $ 偏航舵偏角 $\delta_y$ $ \left [ -20^{\circ},20^{\circ} \right ] $ 俯仰舵偏角 $\delta_z$ $ \left [ -20^{\circ},20^{\circ} \right ] $ 折叠角 $\delta_f$ $\left [ -30^{\circ},155^{\circ} \right ]$ 考虑如下条件: 零舵偏, 即$ {{\delta }_{x}} = {{\delta }_{y}} = {{\delta }_{z}} = {{0}^{\circ }} $; 不考虑变形参数, 零折叠角, 即$ {{\delta }_{f}} = {{0}^{\circ }} $; 零侧滑角, 即$ \beta = {{0}^{\circ }} $. 在上述条件下, 绘制不同马赫数下升阻比随攻角变化的曲线, 如图3所示. 由图可知, 宽速域下, 升阻比呈先增大后减小趋势, 在$ \alpha = {{6}^{\circ }} $处获得最大升阻比. 进一步考虑变形影响, 绘制零舵偏、零侧滑角、${Ma} = 5 $时不同攻角下升阻比随折叠角变化的曲线, 如图4所示. 由图可知, 当折叠角为负时, 随着折叠角逐渐增大并趋近于零, 升阻比呈现增大趋势, 当折叠角为0° 时, 升阻比最大; 当折叠角为正时, 随着折叠角逐渐增大至90°, 升阻比逐渐减小, 折叠角大于90° 时, 升阻比基本维持不变. 图5 ~ 7给出了不同攻角、侧滑角条件下气动力矩系数随折叠角变化的曲线. 对于横侧向力矩系数, 必须考虑侧滑角的影响. 由图5和图6可知, 随着折叠角变化, 滚转力矩系数先增大后减小, 转折折叠角为60°; 而偏航力矩系数变化则更为复杂, 分别在0°、60° 和90° 改变升降趋势. 对于纵向力矩系数, 即俯仰力矩系数, 由式(7) 和图7可知, 其主要受攻角和折叠角影响, 折叠变形的影响特性与升阻比基本保持一致.
图 5 滚转力矩系数随折叠角变化曲线Fig. 5 Curves of rolling moment coefficient varying with folding angle
图 6 偏航力矩系数随折叠角变化曲线Fig. 6 Curves of yawing moment coefficient varying with folding angle
图 7 俯仰力矩系数随折叠角变化曲线Fig. 7 Curves of pitching moment coefficient varying with folding angle通过上述对飞行器气动特性的分析, 可以得到如下主要结论:
1) 对于升阻特性, 当$ \alpha = {{6}^{\circ }} $时, 飞行器升阻比最大; 折叠变形对升阻比影响较大, $ {{\delta }_{f}} = {{0}^{\circ }} $时, 升阻比最大, $ {{\delta }_{f}} = {{155}^{\circ }} $时, 升阻比最小. 因此可以根据需要设计折叠角剖面, 以达到飞行器增程的目的.
2) 对于俯仰力矩系数, 折叠变形对其数值上的影响与升阻比类似; 折叠变形影响$ \frac{\partial {{C}_{mz}}}{\partial \alpha } $, 即影响飞行器静稳定性, 具体来看, 当$ {{\delta }_{f}}<{{30}^{\circ }} $时, 飞行器静不稳定; 当$ {{\delta }_{f}}\ge {{30}^{\circ }} $时, 飞行器静稳定. 显然, 变形会进一步增大姿态控制系统的压力.
3) 对于滚转力矩系数和偏航力矩系数, 折叠变形对其产生不同影响趋势, 同时也与俯仰力矩系数不同. 总的来说, 当$ {{\delta }_{f}}\in \left[ {{0}^{\circ }},{{60}^{\circ }} \right] $, 无论是纵向参数还是横侧向参数, 折叠变形对气动参数的影响成单调递增或递减趋势, 交联耦合较少; 当$ {{\delta }_{f}}\in \left( {{60}^{\circ }},{{90}^{\circ }} \right] $时, 折叠变形对横侧向气动参数影响趋势发生转折, 且远离力矩平衡状态; 当$ {{\delta }_{f}}\in \left( {{90}^{\circ }},{{155}^{\circ }} \right] $时, 折叠变形对气动参数影响不大, 且会导致机体表面气流耦合加剧, 增大模型不确定性和外界干扰. 因此, 在本文中, 仅选用$ {{\delta }_{f}}\in \left[ {{0}^{\circ }},\text{ }{{60}^{\circ }} \right] $ 和 $ {{\delta }_{f}} = {{155}^{\circ }} $时的变形参数, 确保获得优良气动特性的同时减轻姿态控制系统的压力.
4) 侧滑角对滚转力矩系数和偏航力矩系数影响较大, 折叠变形会进一步增大影响趋势, 且增大横纵向气动耦合效应. 本文飞行器采用倾斜转弯(Bank to turn, BTT)方式, 需保证侧滑角为0°, 以提升姿态控制性能.
注 1. 对高超声速变外形飞行器而言, 折叠变形带来的最大影响便是对升阻特性的影响, 变化值最大可达1.4. 在实际中, 可根据不同任务需要调整折叠角, 其带来的机动性和灵活性可有效提高飞行器任务适应性和实现能力. 同时, 机翼折叠造成的升阻特性变化, 可对整个飞行弹道产生巨大影响, 有效提高射程.
注 2. 根据前文结论可知, 变形参数满足${{\delta }_{f}}\in [ {{0}^{\circ }},{{60}^{\circ }} ]$ 或 $ {{\delta }_{f}} = {{155}^{\circ }} $. 这里是指, 只选择$ {{\delta }_{f}}\in [ {{0}^{\circ }},{{60}^{\circ }} ] $和$ {{\delta }_{f}} = {{155}^{\circ }} $时的飞行器构型作为中间形态, 飞行器长时间在该范围内的折叠角剖面下滑翔飞行; ${{\delta }_{f}}\in ( {{60}^{\circ }},{{155}^{\circ }} )$为短期过渡形态, 仅在机翼完全折叠时短暂出现.
1.4 控制问题描述
考虑模型不确定性及外部干扰, 可将式(2)改写为如下形式:
$$ \left\{ \begin{aligned} & {{{\dot{\boldsymbol{X}}}}_{1}} = {{\boldsymbol{X}}_{2}} \\ & {{{\dot{\boldsymbol{X}}}}_{2}} = {{\boldsymbol{f}}_{1}}+{{{\boldsymbol{g}}}_{1}}{\boldsymbol{u}}+{{\boldsymbol{d}}_{1}} \end{aligned} \right. $$ (14) 式中,
$$ \left\{ \begin{aligned} &{{\boldsymbol{X}}_{1}} = {\boldsymbol{\Theta }} \\ & {{\boldsymbol{X}}_{2}} = {\boldsymbol{R}}{\boldsymbol{\omega }} \\ &{{\boldsymbol{f}}_{1}} = {{\dot{\boldsymbol{R}}}{\boldsymbol{\omega}} }-{\boldsymbol{R}}{{\boldsymbol{I}}^{-1}}{{{\boldsymbol{\omega }}}^{\times }}{{\boldsymbol{I}}{\boldsymbol{\omega}} } \\ & {{{\boldsymbol{g}}}_{1}} = {\boldsymbol{R}}{{\boldsymbol{I}}^{-1}} \\ &{{\boldsymbol{d}}_{1}} = {\boldsymbol{R}}{{\boldsymbol{I}}^{-1}}\Delta {{\boldsymbol{d}}_{1}}+\Delta {{\boldsymbol{d}}_{2}} \\ & \Delta {{\boldsymbol{d}}_{2}} = \Delta {{\boldsymbol{f}}_{1}}+\Delta {{{\boldsymbol{g}}}_{1}} +\Delta {{\boldsymbol{f}}_{\delta}}\\ &{\boldsymbol{u}} = {{\boldsymbol{M}}_{t}} = {{{\boldsymbol{f}}}_{\delta }}\left( {{{\boldsymbol{\delta }}}_{E}} \right) \\ &{{{\boldsymbol{\delta }}}_{E}} = {\boldsymbol{f}}_{\delta }^{-1}\left( {\boldsymbol{u}} \right) \end{aligned} \right. $$ (15) 式中, $ {\boldsymbol{\Theta }} = {{[\alpha ,\beta ,\sigma ]}^{\text{T}}} $为姿态角向量; $ {\boldsymbol{u}} $为控制力矩向量, $ {{{\boldsymbol{\delta }}}_{E}} = {{[{{\delta }_{x}},{{\delta }_{y}},{{\delta }_{z}}]}^{\text{T}}} $为等效舵偏角向量, $ {{\boldsymbol{f}}_{\delta}} $为$ {{{\boldsymbol{\delta }}}_{E}} $和$ {{\boldsymbol{M}}_{t}} $的映射函数, 联立式(4) ~ (13)即可唯一确定; $ {{\boldsymbol{d}}_{1}} $为复合总扰动, 由外部干扰力矩$ \Delta {{\boldsymbol{d}}_{1}} $和模型不确定项$ \Delta {{\boldsymbol{d}}_{2}} $构成, $ \Delta {{\boldsymbol{f}}_{1}} $, $ \Delta {{{\boldsymbol{g}}}_{1}} $, $ \Delta{{\boldsymbol{f}}_{\delta}} $表示模型不确定性导致的偏差量; 矩阵$ {\boldsymbol{I}} $, $ {{{\boldsymbol{\omega }}}^{\times}} $和$ {\boldsymbol{R}} $可表示如下:
$$ {{\boldsymbol{I}}}=\left[ \begin{matrix} {{I}_{xx}} & -{{I}_{xy}} & 0 \\ -{{I}_{xy}} & {{I}_{yy}} & 0 \\ 0 & 0 & {{I}_{zz}} \\ \end{matrix} \right] $$ (16) $$ {{\boldsymbol{\omega }}^{\times }}=\left[ \begin{matrix} 0 & -{{\omega }_{z}} & {{\omega }_{y}} \\ {{\omega }_{z}} & 0 & -{{\omega }_{x}} \\ -{{\omega }_{y}} & {{\omega }_{x}} & 0 \end{matrix} \right] $$ (17) $$ {{\boldsymbol{R}}}=\left[ \begin{matrix} -\tan \beta \cos \alpha & \tan \beta \sin \alpha & 1 \\ \sin \alpha & \cos \alpha & 0 \\ \sec \beta \cos \alpha & -\sec \beta \sin \alpha & 0 \end{matrix} \right] $$ (18) 本文所研究的姿态控制问题可总结如下:
1) 控制目标. 在飞行器存在总扰动$ {{\boldsymbol{d}}_{1}} $的情况下, 设计干扰观测器以实现对$ {{\boldsymbol{d}}_{1}} $的固定时间精确估计, 并基于估计结果设计控制力矩$ {\boldsymbol{u}} $, 继而得到舵偏角控制量$ {{{\boldsymbol{\delta }}}_{E}} $, 使姿态角${\boldsymbol{\Theta }}$跟踪期望姿态角${{{\boldsymbol{\Theta }}}_{d}}$, 满足预先设定的瞬态性能和稳态性能并能适应不同的变形条件.
2) 扰动来源及形式. 由式(15)可知, 待估计总扰动为$ {{\boldsymbol{d}}_{1}} $, 一般由外部干扰和模型不确定性构成. 根据实际经验, 飞行器所受外部干扰力矩$ \Delta {{\boldsymbol{d}}_{1}} $主要受气流影响. 此外, 飞行器折叠翼大尺度变形会对惯量参数和等效机体参数产生较大影响. 对转动惯量而言, 一般难以准确获得变形过程中惯量参数变化的实时特性, 但结合第1.3节飞行器气动特性分析结论和式(8)可知, 惯量参数主要与$ \sin \delta_f $或$ \cos \delta_f $有关, 故可采用式(19)对飞行器转动惯量进行近似计算.
$$ \left\{\begin{aligned} &I_{xx} = I_{xx,0}+\frac{(I_{xx,1}-I_{xx,0})(\cos\delta_{f,1}-\cos\delta_{f})}{\cos\delta_{f,1}-\cos\delta_{f,0}}\\ & I_{yy} = I_{yy,0}+\frac{(I_{yy,1}-I_{yy,0})(\sin\delta_{f,1}-\sin\delta_{f})}{\sin\delta_{f,1}-\sin\delta_{f,0}}\\ & I_{zz} = I_{zz,0}+\frac{(I_{zz,1}-I_{zz,0})(\cos\delta_{f,1}-\cos\delta_{f})}{\cos\delta_{f,1}-\cos\delta_{f,0}}\\ &I_{xy} = I_{xy,0}+\frac{(I_{xy,1}-I_{xy,0})(\sin\delta_{f,1}-\sin\delta_{f})}{\sin\delta_{f,1}-\sin\delta_{f,0}} \end{aligned}\right. $$ (19) 式中, $ I_{xx,0} $, $ I_{yy,0} $, $ I_{zz,0} $, $ I_{xy,0} $分别为转动惯量的下界, 对应$ \delta_{f,1} = 155^{\circ} $; $ I_{xx,1} $, $ I_{yy,1} $, $ I_{zz,1} $, $ I_{xy,1} $分别为转动惯量的上界, 对应$ \delta_{f,0} = 0^{\circ} $. 因此, 需要考虑上述惯量参数的估计误差$ \Delta {\boldsymbol{I}} $, 此为偏差量$ \Delta {{\boldsymbol{f}}_{1}} $, $ \Delta {{{\boldsymbol{g}}}_{1}} $的主要来源. 机翼折叠变形还会对气动力系数、气动力矩系数、质心和压心位置以及等效参考长度等参数造成不确定影响, 进而引入偏差量$ \Delta{{\boldsymbol{f}}_{\delta}} $. 综上可知, 总扰动$ {{\boldsymbol{d}}_{1}} $包含了复杂多源的不确定性, 且同样具有快时变、强耦合和强非线性的特点, 有必要对其进行快速精确估计, 从而增强姿态控制性能. 在本文中, 所考虑总扰动$ {{\boldsymbol{d}}_{1}} $始终满足假设4.
3) 参考指令. 飞行器参考指令主要包括期望姿态角${{{\boldsymbol{\Theta }}}_{d}}$和期望折叠角$ \delta_f $. ${{\boldsymbol{\Theta}}_{d}} = {{[\alpha_d ,\beta_d ,\sigma_d ]}^{\text{T}}}$一般由制导回路给出. 本文采取开环独立变形控制策略, 即折叠翼的变形方案由飞行任务的先验内容提前确定, 此时的控制方案类似于开环控制, 飞行控制系统的作用是实现变外形飞行器的变形与飞行协调控制. 所述折叠翼的变形方案即为折叠角剖面的预先设计. 在本文中, 参考指令始终满足假设5.
4) 执行机构驱动策略. 实际中, 飞行器跨大空域、宽速域飞行, 全飞行包线内需采用反作用控制系统(Reaction control system, RCS)和气动舵等多种执行机构来进行复合控制. 考虑到本文仅研究30 ~ 40 km高度无动力滑翔段的飞行器姿态控制问题, 故仅以气动舵为执行机构. 根据图1可知, 飞行器采用4片“十”字形配置尾舵作为控制舵, 尾舵实际舵偏角向量$ {{\boldsymbol{\delta}}_{R}} = {{[\delta_1, \delta_2, \delta_3, \delta_4]}^{\text{T}}} $, 且$ {{\boldsymbol{\delta}}_{E}} $ = $ {{\boldsymbol{M}}_{\delta}}{{\boldsymbol{\delta}}_{R}} $, $ {{\boldsymbol{M}}_{\delta}} $为舵面控制分配矩阵, 可表示如下:
$$ {{\boldsymbol{M}}_{\delta}}=\left[ \begin{array}{*{20}{c}} \dfrac{1}{4} & \dfrac{1}{4} & -\dfrac{1}{4} & -\dfrac{1}{4} \\ \dfrac{1}{2} & 0 &\dfrac{1}{2} & 0 \\ 0 & \dfrac{1}{2} & 0 & \dfrac{1}{2} \end{array} \right] $$ (20) 假设 4[2, 9]. 总扰动$ {{d}_{1, i}} $有界、连续、一阶可导, 且满足: $ \left| {{d}_{1,i}} \right|\le {{L}_{d1}} $, $ | {{\dot{d}}_{1,i}} |\le {{L}_{d2}} $, 其中, $ {{L}_{d1}} $, $ {{L}_{d2}} $为已知常数.
假设 5[22]. 参考指令$ \alpha_d $, $ \beta_d $, $ \sigma_d $, $ \delta_f $均有界、连续、一阶可导, 其一阶导数也有界、连续.
针对上述控制问题, 本文提出了一种基于固定时间预设性能的姿态控制方法, 具体控制方案如图8所示.
图 8 飞行器固定时间预设性能控制方案框图Fig. 8 Flowchart of fixed-time prescribed performance control for HMV注 3. 在本文中, 考虑到飞行器模型中$ {{\boldsymbol{f}}_{\delta}} $函数的复杂性, 采用式(15)表示形式可有效避免$ {{\boldsymbol{f}}_{\delta}} $函数冗长的显式表达, 有利于简化后续控制器的设计流程. 同时, 基于式(4) ~ (13)以及式(20)可确定$ {{\boldsymbol{M}}_{t}} $, $ {{\boldsymbol{\delta}}_{R}} $, 和$ {{\boldsymbol{\delta}}_{E}} $三者之间唯一的对应关系. 因此, 以控制力矩作为待设计控制量并解算得到舵偏角具有充分可行性.
2. 基于预设性能的固定时间姿态控制方法
2.1 预备知识
为便于后续控制方法设计, 下面将介绍本文所用到的定义及引理.
定义 1. 1) $ {{\bf{R}}^{n}} $表示$ n $维欧氏空间, $ {{\bf{R}}_{+}} $表示正实数; 2) 对于向量$ {\boldsymbol{a}} = {{\left[ {{a}_{1}}, {{a}_{2}}, {{a}_{3}} \right]}^{\text{T}}} $和$ {\boldsymbol{b}}\in {{\bf{R}}^{3}} $, $ {\boldsymbol{a}}\times {\boldsymbol{b}} $表示两向量叉乘, $ {\boldsymbol{a}}\times {\boldsymbol{b}} = {\boldsymbol{S}}({\boldsymbol{a}}){\boldsymbol{b}} $, 其中,
$$ {{\boldsymbol{S}}}({{\boldsymbol{a}}})=\left[ \begin{matrix} 0 & -{{a}_{3}} & {{a}_{2}} \\ {{a}_{3}} & 0 & -{{a}_{1}} \\ -{{a}_{2}} & {{a}_{1}} & 0 \end{matrix} \right] $$ 3) 对任意向量${{\boldsymbol{a}}_{*}}\; = \;{{\left[ {{a}_{*,1}},{{a}_{*,2}},\cdots ,{{a}_{*,n}} \right]}^{\text{T}}}$, ${{a}_{*,i}}\text{ }(i = 1,2,\cdots ,n)$为$ {{\boldsymbol{a}}_{*}} $的分量形式, 当$ n = 3 $时, 除非特别说明, 一般省略其后缀$ (i = 1,2,3) $, $ \left\| {{\boldsymbol{a}}_{*}} \right\| $为向量$ {{\boldsymbol{a}}_{*}} $的二范数; 4) 对$ x\in {\bf{R}} $, $ p\in {\bf{R}} $, 令${{\left\lfloor x \right\rfloor }^{p}} = \text{sgn} (x){{\left| x \right|}^{p}}$, $ \text{sgn}(x) $为符号函数.
定义 2[23]. 令$ {\cal{D}} $为一包含原点的开区域, 积分障碍Lyapunov函数(Integral barrier Lyapunov function, IBLF) $ V(x) $为定义在$ {\cal{D}} $上关于系统$ \dot{x} = f(x) $的标量函数, 它具有以下特性: 1) 光滑、正定; 2) $ {\cal{D}} $上的每个点的一阶连续偏导数存在; 3) 当$ x $趋向$ {\cal{D}} $的边界时, $ V(x)\to \infty $; 4) 当$ x(0)\in {\cal{D}} $时, $ \forall t\ge 0,\text{ }V(x(t))\le b $, 其中, $ b $为某个正常数.
引理 1[24]. 对于初始条件有界的系统, 若存在一个 $ {{\cal{C}}^{1}} $连续且正定的Lyapunov函数$ V(x) $, 满足$ \left\| {{\pi }_{1}} \right\|\le V(x)\le \left\| {{\pi }_{2}} \right\| $, 若$ \dot{V}(x) \le {{c}_{1}}V(x)+{{c}_{2}} ,$ 则系统的解$ x(t) $一致有界, 其中, $ {{\pi }_{1}} $, $ {{\pi }_{2}} $为$ {{\kappa }_{\infty }} $类函数, $ {{c}_{1}} $, $ {{c}_{2}} $为正常数.
引理 2[25]. 定义开集$ {{\bf{Q}}_{1}} $及区间$ {{\bf{N}}_{1}} $:
$$ \left\{ \begin{aligned} & {{\bf{Q}}_{1}}: = \left\{ {{q}_{1}}\in {\bf{R}}:-{{k}_{b}}<{{q}_{1}}<{{k}_{b}} \right\}\subset {\bf{R}},\text{ }{{k}_{b}}\in {{\bf{R}}_+} \\ & {{\bf{N}}_{1}}: = {{\bf{R}}^{l}}\times {{\bf{Q}}_{1}}\subset {{\bf{R}}^{l+1}},\text{ }l\in {{\bf{R}}_+} \end{aligned} \right. $$ 考虑系统
$$ {{\dot{\boldsymbol{Q}}}_{0}} = {{h}_{1}}(t,{{\boldsymbol{Q}}_{0}}) $$ (21) 式中, $ {{\boldsymbol{Q}}_{0}}: = {{\left[ {{q}_{1}},{{q}_{2}} \right]}^{\text{T}}} $$ \in {\bf{N}_1} $与$ {{h}_{1}}: = {{\bf{R}}_{+}}\times {{\bf{N}}_{1}} $$ \to $$ {{\bf{R}}^{l+1}} $分段连续于$ t $, 局部Lipschitz于$ {{q}_{1}} $, 且一致于$ t $和$ {{\bf{R}}_{+}}\times{{\bf{N}}_{1}} $. 假设存在连续可微的正定函数$ {{V}_{L1}}: = {{\bf{Q}}_{1}} \to {{\bf{R}}_{+}} $ 和 $ {{V}_{L2}}: = {{\bf{R}}^{l}}\to {{\bf{R}}_{+}} $, 满足:
$$ \left\{ \begin{aligned} & \underset{|{{q}_{1}}|\to {{k}_{b}}}{\mathop{\lim }}\,{{V}_{L1}}({{q}_{1}}) = \infty \\ & {{\pi }_{3}}(\left\| {{q}_{2}} \right\|)\le {{V}_{L2}}({{q}_{2}})\le {{\pi }_{4}}(\left\| {{q}_{2}} \right\|) \end{aligned} \right. $$ (22) 式中, $ {{\pi }_{3}} $和$ {{\pi }_{4}} $ 均为$ {{\kappa }_{\infty }} $类函数. 令$V({{\boldsymbol{Q}}_{0}}) = {{V}_{L1}}({{q}_{1}})\;+ {{V}_{L2}}({{q}_{2}})$, $ {{q}_{1}}(0)\in {{\bf{Q}}_{1}} $, 若满足:
$$ \dot{V} = \frac{\partial V}{\partial {{\boldsymbol{Q}}_{0}}}{{h}_{1}}\le -{{c}_{3}}V+{{c}_{4}} $$ (23) 则$ \text{ }\forall t\in \left[ 0,\infty \right) $, $ {{q}_{1}}(t)\in {{\bf{Q}}_{1}} $, 其中, $ {{c}_{3}} $, $ {{c}_{4}} $为正常数.
2.2 固定时间干扰观测器设计
基于文献[26], 本节设计了一种固定时间干扰观测器, 可实现对总扰动$ {{\boldsymbol{d}}_{1}} $的固定时间精确估计. 观测器形式如下:
$$ \left\{ \begin{aligned} {{{\dot{z}}}_{1,i}}& = {{z}_{0,i}}+\frac{{{k}_{z1,i}}}{{{\varepsilon }_{0,i}}}{{\varphi }_{1}}\left( \frac{{{X}_{2,i}}-{{z}_{1,i}}}{{{\varepsilon }_{0,i}}} \right)+{{z}_{2,i}} \\ {{{\dot{z}}}_{2,i}}& = \frac{{{k}_{z2,i}}}{{{\varepsilon }_{0,i}}}{{\varphi }_{2}}\left( \frac{{{X}_{2,i}}-{{z}_{1,i}}}{{{\varepsilon }_{0,i}}} \right) \end{aligned} \right. $$ (24) 式中, $ {{\boldsymbol{z}}_{0}} = {{\boldsymbol{f}}_{1}}+{{{\boldsymbol{g}}}_{1}}{\boldsymbol{u}} $, $ {{z}_{1,i}} $, $ {{z}_{2,i}} $为观测器状态量; 修正项$ {{\varphi}_{1}}\left( x \right) $, $ {{\varphi}_{2}}\left( x \right) $的形式如下:
$$ \left\{ \begin{aligned} & {{\varphi }_{1}}\left( x \right) = {{\left\lfloor x \right\rfloor }^{{{\gamma }_{1}}}}+{{\left\lfloor x \right\rfloor }^{{{\gamma }_{2}}}} \\ & {{\varphi }_{2}}\left( x \right) = {{\left\lfloor x \right\rfloor }^{2{{\gamma }_{1}}-1}}+{{\left\lfloor x \right\rfloor }^{2{{\gamma }_{2}}-1}} \end{aligned} \right. $$ (25) $ {{k}_{z1,i}} $, $ {{k}_{z2,i}} $为观测器增益, 满足
$$ \left\{ {{k}_{z1,i}},{{k}_{z2,i}}\in {{\bf{R}}_+}\left| {{k}_{z1,i}}\ge 2\sqrt{{{k}_{z2,i}}} \right. \right\} $$ (26) $ {{\varepsilon }_{0,i}}\in \left( 0,1 \right) $ 为放大因子; $ {{\gamma }_{1}}\in \left( 0.5,1 \right) $, $ {{\gamma }_{2}}\in \left( 1,1.5 \right) $. 定义观测误差:
$$ \left\{ \begin{aligned} & {{e}_{z1,i}} = {{X}_{2,i}}-{{z}_{1,i}} \\ & {{e}_{z2,i}} = {{d}_{1,i}}-{{z}_{2,i}} \end{aligned} \right. $$ (27) 令$ {{z}_{\eta j,i}} = {{{e}_{zj,i}}\left( {{\varepsilon }_{0,i}}t \right)}/{\varepsilon _{0,i}^{2-j}} $, $ j = 1,2 $, 可得
$$ \left\{ \begin{aligned} {{{\dot{z}}}_{\eta 1,i}} = {}&-{{k}_{z1,i}}\left( {{\left\lfloor {{z}_{\eta 1,i}} \right\rfloor }^{{{\gamma }_{1}}}}+{{\left\lfloor {{z}_{\eta 1,i}} \right\rfloor }^{{{\gamma }_{2}}}} \right)+{{z}_{\eta 2,i}} \\ {{{\dot{z}}}_{\eta 2,i}} = {}&-{{k}_{z2,i}}\left( {{\left\lfloor {{z}_{\eta 1,i}} \right\rfloor }^{2{{\gamma }_{1}}-1}}+{{\left\lfloor {{z}_{\eta 1,i}} \right\rfloor }^{2{{\gamma }_{2}}-1}} \right)+\\ {}&{{\varepsilon }_{0,i}}{{{\dot{d}}}_{1,i}} \end{aligned} \right. $$ (28) 定义$ {{\boldsymbol{z}}_{\eta ,i}} = {{\left[ {{z}_{\eta 1,i}},{{z}_{\eta 2,i}} \right]}^{\text{T}}} $, $ {{d}_{2,i}} = {{\varepsilon }_{0,i}}{{\dot{d}}_{1,i}} $.
根据文献[26]的定理1可知, 对于系统(28), 误差向量$ {{\boldsymbol{z}}_{\eta ,i}} $将在固定时间内收敛至原点的某个邻域内, 收敛域可表示如下:
$$ \left\{ \begin{aligned} & {{{\boldsymbol{E}}}_{1,i}} = \left\{ \left. {{\boldsymbol{z}}_{\eta ,i}}\in {{\bf{R}}^{2}} \right|{{V}_{z}}\left( {{\boldsymbol{z}}_{\eta ,i}} \right)\le {{V}_{z1,i}}\right\},\text{ }{{d}_{3,i}}>1 \\ & {{{\boldsymbol{E}}}_{2,i}} = \left\{ \left. {{\boldsymbol{z}}_{\eta ,i}}\in {{\bf{R}}^{2}} \right|{{V}_{z}}\left( {{\boldsymbol{z}}_{\eta ,i}} \right)\le {{V}_{z2,i}}\right\},\text{ }{{d}_{3,i}}\le 1 \end{aligned} \right. $$ (29) $$ \left\{ \begin{aligned} & {{V}_{z1,i}}: = {{\left( {{d}_{3,i}} \right)}^{\frac{2{{\gamma }_{1}}}{2{{\gamma }_{1}}-1}}} \\ & {{V}_{z2,i}}: = {{\left( {{d}_{3,i}} \right)}^{\frac{2{{\gamma }_{2}}}{2{{\gamma }_{2}}-1}}} \end{aligned} \right. $$ (30) 收敛时间上界可表示如下:
$$ \left\{ \begin{aligned} & {{T}_{z1,i}} = \frac{2{{\gamma }_{2}}}{{{c}_{01}}\left( {{\gamma }_{2}}-1 \right)}{{\left( {{d}_{3,i}} \right)}^{\frac{-{{\gamma }_{1}}\left( {{\gamma }_{2}}-1 \right)}{{{\gamma }_{2}}\left( 2{{\gamma }_{1}}-1 \right)}}},& \text{ }{{d}_{3,i}}>1 \\ & {{T}_{z2,i}} = \frac{1}{{{c}_{01}}}\left( \frac{2{{\gamma }_{1}}}{\left( 1-{{\gamma }_{1}} \right)}+\frac{2{{\gamma }_{2}}}{\left( {{\gamma }_{2}}-1 \right)} \right),& \text{ }{{d}_{3,i}}\le 1 \end{aligned} \right. $$ (31) 式中, $ {{d}_{3,i}} = {2{{c}_{02}}{{d}_{2,i}}}/{{{c}_{01}}} $, $ {{c}_{01}},{{c}_{02}}\in {{\bf{R}}_{+}} $为常数.
综上所述, 当总扰动满足假设3条件时, 采用式(24)所示干扰观测器, 可保证$ {{{z}}_{2,i}} $在固定时间$ {{T}_{z1,i}} $, $ {{T}_{z2,i}} $内实现对$ {{{d}}_{1,i}} $的精确估计.
注 4. 对于式(28)所示观测器系统, 当参数$ {{c}_{1}} $, $ {{c}_{2}} $确定后, 若$ {{\varepsilon }_{0,i}}\to 0 $, 则$ {{V}_{z1,i}},{{V}_{z2,i}}\to 0 $, 即干扰观测器系统可在理论上实现精确收敛. 与传统的扰动估计方法如扩张状态观测器[1]、干扰观测器[13]不同, 本部分采用的干扰观测器可以实现固定时间精确收敛, 适用于一些需要快速响应的实时任务.
2.3 固定时间预设性能函数设计
为定量描述姿态跟踪性能, 本节基于固定时间控制理论建立了一个新型固定时间预设性能函数(Prescribed performance function, PPF). 该性能函数的形式及性质由定理1给出.
定理 1. 考虑性能函数$ \rho (t) $, 其一阶导数为
$$ \dot{\rho }(t) = \left\{ \begin{aligned} & -{{\alpha }_{0}}\left( {{\alpha }_{01}}{{\rho }_{e}}{{(t)}^{{{\sigma }_{1}}}}+{{\alpha }_{02}}{{\rho }_{e}}{{(t)}^{{{\sigma }_{2}}}} \right),&\text{ }{{\rho }_{e}}(t)>0 \\ & 0,&\text{ }{{\rho }_{e}}(t) = 0 \end{aligned} \right. $$ (32) 式中,
$$ \left\{ \begin{aligned} & {{\sigma }_{1}} = {{\left( \frac{{{m}_{1}}}{{{m}_{2}}}\; \right)}^{\text{sgn} \left( 1-\left| {{\rho }_{e}}(t) \right| \right)}} \\ &{{\sigma }_{2}} = {{\left( \frac{{{n}_{1}}}{{{n}_{2}}}\; \right)}^{\text{sgn} \left( 1-\left| {{\rho }_{e}}(t) \right| \right)}} \\ & {{\rho }_{e}}(t) = \rho (t)-{{\rho }_{\infty }} \end{aligned} \right. $$ (33) 式中, $ {{\alpha }_{0}} $, $ {{\alpha }_{01}} $, $ {{\alpha }_{02}}\in {{\bf{R}}_{+}} $均为常数, 且$ {{\alpha }_{0}}{{\alpha }_{01}} = {{\alpha }_{1}} $, $ {{\alpha }_{0}}{{\alpha }_{02}} = {{\alpha }_{2}} $; $ {{m}_{1}} $, $ {{m}_{2}} $, $ {{n}_{1}} $, $ {{n}_{2}}\in {{\bf{R}}_{+}} $为整奇数, 且满足$ {{m}_{1}}<{{m}_{2}} $, $ {{n}_{1}}<{{n}_{2}} $; $ {{\rho }_{0}} $和$ {{\rho }_{\infty }} $分别为性能函数初值和终值, 即$ \rho (t = 0) = {{\rho }_{0}} $, $ \rho (t\to \infty ) = {{\rho }_{\infty }} $, 且$ {{\rho }_{0}}> {{\rho }_{\infty }}>0 $. 此时, $ \rho (t) $可在时间$ {{T}_{0}} $内收敛至$ {{\rho }_{\infty }} $, 并在$ t>{{T}_{0}} $后保持不变, 且$ {{T}_{0}} $满足:
$$ \begin{split}{{T}_{0}}<{{T}_{0,\max }} =\;& \min \left\{ {{T}_{01,m}},{{T}_{01,n}} \right\}+\\ &\min \left\{ {{T}_{02,m}},{{T}_{02,n}} \right\} \end{split}$$ (34) 证明. 建立如下Lyapunov函数:
$$ {{V}_{0}} = \frac{1}{2}{{\left( \rho (t)-{{\rho }_{\infty }} \right)}^{2}} $$ (35) 对$ {{V}_{0}} $求导可得
$$ \begin{split} {{{\dot{V}}}_{0}} =\; & \left( \rho (t)-{{\rho }_{\infty }} \right)\dot{\rho }(t) = \\ & -{{\alpha }_{1}}{{\left( 2{{V}_{0}} \right)}^{\frac{1+{{\sigma }_{1}}}{2}}}-{{\alpha }_{2}}{{\left( 2{{V}_{0}} \right)}^{\frac{1+{{\sigma }_{2}}}{2}}} \end{split} $$ (36) 由式(36)可知, $ {{\dot{V}}_{0}} $负定, 因此$ \rho (t) $是稳定的. 基于文献[27], 为进一步获得$ \rho (t) $收敛时间的上界, 可将其收敛过程划分为两个阶段. 因此, 式(32)可改写为
$$ \dot{\rho }(t) = \left\{ \begin{aligned} &-{{\alpha }_{1}}{{\rho }_{e}}{(t)}^{\frac{m_2}{m_1}}-{{\alpha }_{2}}{{\rho }_{e}}{(t)}^{\frac{n_2}{n_1}},\text{ }{{\rho }_{e}}(t)\ge 1 \\ &-{{\alpha }_{1}}{{\rho }_{e}}{(t)}^{\frac{m_1}{m_2}}-{{\alpha }_{2}}{{\rho }_{e}}{(t)}^{\frac{n_1}{n_2}},\text{ }0\le {{\rho }_{e}}(t)<1 \end{aligned} \right. $$ (37) 1) 当${{\rho }_{e}}(t) \ge 1$时, 定义辅助变量${{y}_{1}} = {{\rho }_{e}}{{(t)}^{1-{{{m}_{2}}}/{{{m}_{1}}}}}$, 考虑到$ 1-{{{m}_{2}}}/{{{m}_{1}}}<0 $, 因此 $ {{y}_{1}}\in \left( 0,1 \right] $, $ {{\rho}_{e}}(t) = {{y}_{1}}^{\frac{{{m}_{1}}}{{{m}_{1}}-{{m}_{2}}}} $, 则可得
$$ {{\dot{y}}_{1}} = -\frac{{{m}_{1}}-{{m}_{2}}}{{{m}_{1}}}\left( {{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}^{\frac{{{n}_{2}}{{m}_{1}}-{{n}_{1}}{{m}_{2}}}{{{n}_{1}}\left( {{m}_{1}}-{{m}_{2}} \right)}} \right) $$ (38) 令$ {{\beta }_{1}} = \frac{{{n}_{2}}{{m}_{1}}-{{n}_{1}}{{m}_{2}}}{{{n}_{1}}\left( {{m}_{1}}-{{m}_{2}} \right)} $, 则$ {{\beta }_{1}}-1 = \frac{{{m}_{1}}\left( {{n}_{2}}-{{n}_{1}} \right)}{{{n}_{1}}\left( {{m}_{1}}-{{m}_{2}} \right)}<0 $, 即$ {{\beta }_{1}}<1 $. 式(38) 可进一步改写为
$$ \begin{split} &\frac{{{m}_{1}}}{{{m}_{2}}-{{m}_{1}}}\frac{\text{d}{{y}_{1}}}{\text{d}t} = {{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}^{{{\beta }_{1}}}\text{ } \Leftrightarrow \\ &\qquad \frac{1}{{{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}^{{{\beta }_{1}}}}\text{d}{{y}_{1}} = \frac{{{m}_{2}}-{{m}_{1}}}{{{m}_{1}}}\text{d}t \end{split} $$ (39) 然后, 对式(39)两端分别进行积分, 可得
$$ \varphi ({{y}_{1}}) = \int_{0}^{{{y}_{1}}}{\frac{1}{{{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}^{{{\beta }_{1}}}}\text{d}{{y}_{1}}} = \frac{{{m}_{2}}-{{m}_{1}}}{{{m}_{1}}}t $$ (40) 式中, $ \varphi ({{y}_{1}}) $在$ {{y}_{1}}\in \left( 0,1 \right] $上单调递增, 当且仅当$ t = {{T}_{01}} $时, $ {{y}_{1}} = 1 $, 因此 $ {{T}_{01}} $可由式(41)得到.
$$ \begin{split} {{T}_{01}} =\; &\frac{{{m}_{1}}}{{{m}_{2}}-{{m}_{1}}}\int_{0}^{1}{\frac{1}{{{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}^{{{\beta }_{1}}}}\text{d}{{y}_{1}}} <\\ & \frac{{{m}_{1}}}{{{m}_{2}}-{{m}_{1}}}\int_{0}^{1}{\frac{1}{{{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{1}}}\text{d}{{y}_{1}}} = \\ & \frac{{{m}_{1}}}{{{m}_{2}}-{{m}_{1}}}\frac{1}{{{\alpha }_{2}}}\ln \left( 1+\frac{{{\alpha }_{2}}}{{{\alpha }_{1}}} \right): = {{T}_{01,m}} \end{split} $$ (41) 同理, 若定义辅助变量$ {{y}_{1}} = {{\rho }_{e}}{{(t)}^{1-{{{n}_{2}}}/{{{n}_{1}}}}} $, 则
$$ {{T}_{01}}<\frac{{{n}_{1}}}{{{n}_{2}}-{{n}_{1}}}\frac{1}{{{\alpha }_{1}}}\ln \left( 1+\frac{{{\alpha }_{1}}}{{{\alpha }_{2}}} \right): = {{T}_{01,n}} $$ (42) 2) 当$ 0\;\le \;{{\rho }_{e}}(t)\;<\;1 $时, 定义辅助变量 ${{y}_{2}}\; = {{\rho }_{e}}{{(t)}^{1-{{{m}_{1}}}/{{{m}_{2}}}}}$, 考虑到$ 0<1-{{{m}_{1}}}/{{{m}_{2}}}<1 $, 因此, $ {{y}_{2}}\in \left[ 0,1 \right) $, $ {{\rho}_{e}}(t) = {{y}_{2}}^{\frac{{{m}_{2}}}{{{m}_{2}}-{{m}_{1}}}} $. 与式(38) ~ (40)推导过程类似, 可得如下关系:
$$ \varphi ({{y}_{2}}) = \int_{1}^{{{y}_{2}}}{\frac{1}{{{\alpha }_{1}}+{{\alpha }_{2}}{{y}_{2}}^{{{\beta }_{2}}}}\text{d}{{y}_{2}}} = \frac{{{m}_{1}}-{{m}_{2}}}{{{m}_{2}}}t $$ (43) 式中, $ {{\beta }_{2}} = \frac{{{n}_{1}}{{m}_{2}}-{{n}_{2}}{{m}_{1}}}{{{n}_{2}}\left( {{m}_{2}}-{{m}_{1}} \right)}<1 $, $ \varphi ({{y}_{2}}) $在定义域$ {{y}_{2}}\in [ 0, 1 ) $上单调递减, 当且仅当 $ t = {{T}_{01}}+{{T}_{02}} $时, $ {{y}_{2}} = 0 $, 因此$ {{T}_{02}} $可由式(44)得到.
$$ \begin{split} {{T}_{02}} =\; &\frac{{{m}_{2}}}{{{m}_{2}}-{{m}_{1}}}\int_{1}^{0}{\frac{1}{-{{\alpha }_{1}}-{{\alpha }_{2}}{{y}_{2}}^{{{\beta }_{2}}}}\text{d}{{y}_{2}}}< \\ & \frac{{{m}_{2}}}{{{m}_{2}}-{{m}_{1}}}\int_{1}^{0}{\frac{1}{-{{\alpha }_{1}}-{{\alpha }_{2}}{{y}_{2}}}\text{d}{{y}_{2}}} = \\ & \frac{{{m}_{2}}}{{{m}_{2}}-{{m}_{1}}}\frac{1}{{{\alpha }_{2}}}\ln \left( 1+\frac{{{\alpha }_{2}}}{{{\alpha }_{1}}} \right): = {{T}_{02,m}} \end{split} $$ (44) 同理, 若定义辅助变量$ {{y}_{2}} = {{\rho }_{e}}{{(t)}^{1-{{{n}_{1}}}/{{{n}_{2}}}}} $, 则
$$ {{T}_{02}}<\frac{{{n}_{2}}}{{{n}_{2}}-{{n}_{1}}}\frac{1}{{{\alpha }_{1}}}\ln \left( 1+\frac{{{\alpha }_{1}}}{{{\alpha }_{2}}} \right): = {{T}_{02,n}} $$ (45) 因此, 结合式(41)、式(42)、式(44)和式(45), 可得
$$ \begin{split} {{T}_{0}} = \;& {{T}_{01}}+{{T}_{02}}< \\ {}& \left( \min \left\{ {{T}_{01,m}},{{T}_{01,n}} \right\}+\right. \\ {}& \left. \min \left\{ {{T}_{02,m}},{{T}_{02,n}} \right\} \right): = {{T}_{0,\max }} \end{split} $$ (46) 综上, 当$ t\le {{T}_{0}} $时, $\mathop {\lim }\nolimits_{t \le {T_0}} \,\rho (t) = {{\rho }_{\infty }}$; 当$ t>{{T}_{0}} $时, $ \dot{\rho }(t) = 0 $, 即$ \rho (t) $ 保持不变.
□ 注 5. 通过本文所提出的预设性能函数, 可实现在预期固定时间内完成控制任务. 图9给出了不同预期收敛时间上界及稳态误差边界的性能函数变化曲线, 由图可知, 可根据实际控制任务, 预先设定性能函数的收敛特性, 以定量约束姿态跟踪的瞬态性能和稳态性能.
2.4 预设性能姿态控制器设计
在本节中, 为定量描述姿态跟踪性能, 定义如下预设性能约束:
$$ -{{\delta }_{1,i}}{{\rho }_{1,i}}(t)<{{e}_{1,i}}<{{\delta }_{2,i}}{{\rho }_{1,i}}(t) $$ (47) 式中, $ {{\rho}_{1,i}}(t) $为前文所定义性能函数; $ {{\delta }_{1,i}} $, $ {{\delta }_{2,i}}\in {{\bf{R}}_{+}} $为常数; 为简便描述, 定义$ {{\rho }_{i}}(t) = {{\rho }_{i}} $, 即省略下标$ (t) $; $ {{e}_{1,i}} $为定义的姿态角跟踪误差:
$$ {{e}_{1,i}} = {{X}_{1,i}}-{{X}_{1d,i}} $$ (48) 接下来基于动态面控制技术设计姿态跟踪控制器.
步骤 1. 基于定义2, 定义积分障碍Lyapunov函数如下:
$$ \left\{ \begin{aligned} {}&{{V}_{1}} = \sum\limits_{i = 1}^{3}{{{V}_{B,i}}({{e}_{1,i}},{{\rho }_{1,i}})} \\ {}&{{V}_{B,i}}({{e}_{1,i}},{{\rho }_{1,i}}) = (1-{{\ell }_{1}})\int_{0}^{{{e}_{1,i}}}{\frac{s\delta _{1,i}^{2}\rho _{1,i}^{2}}{\delta _{1,i}^{2}\rho _{1,i}^{2}-{{s}^{2}}}{\rm{d}}s}\;+ \\ {}&\qquad{{\ell }_{1}}\int_{0}^{{{e}_{1,i}}}{\frac{s\delta _{2,i}^{2}\rho _{1,i}^{2}}{\delta _{2,i}^{2}\rho _{1,i}^{2}-{{s}^{2}}}{\rm{d}}s} \end{aligned} \right. $$ (49) 式中,
$$ {{\ell }_{1}} = \left\{ \begin{aligned} & 1,\text{ }&{{e}_{1,i}}\ge 0 \\ & 0,\text{ }&{{e}_{1,i}}<0 \end{aligned} \right. $$ (50) 对$ {{V}_{1}} $求导可得
$$ \begin{split} {{{\dot{V}}}_{1}} =\; &\sum\limits_{i = 1}^{3}{\left( {{\eta }_{1,i}}+{{\lambda }_{1,i}}{{{\dot{e}}}_{1,i}} \right)} = \\ &\sum\limits_{i = 1}^{3}{\left( {{\eta }_{1,i}}+{{\lambda }_{1,i}}\left( {{X}_{2,i}}-{{{\dot{X}}}_{1d,i}} \right) \right)} \end{split} $$ (51) 式中,
$$ \left\{ \begin{aligned} {{\lambda }_{1,i}} = {}&(1-{{\ell }_{1}})\frac{{{e}_{1,i}}\delta _{1,i}^{2}\rho _{1,i}^{2}}{\delta _{1,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2}}+{{\ell }_{1}}\frac{{{e}_{1,i}}\delta _{2,i}^{2}\rho _{1,i}^{2}}{\delta _{2,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2}} \\ {{\eta }_{1,i}} = {}&(-1+{{\ell }_{1}})\int_{0}^{{{e}_{1,i}}}{\frac{2{{s}^{3}}\delta _{1,i}^{2}{{\rho }_{1,i}}{{{\dot{\rho }}}_{1,i}}}{{{\left( \delta _{1,i}^{2}\rho _{1,i}^{2}-{{s}^{2}} \right)}^{2}}}\text{d}s}\;- \\ {}& {{\ell }_{1}}\int_{0}^{{{e}_{1,i}}}{\frac{2{{s}^{3}}\delta _{2,i}^{2}{{\rho }_{1,i}}{{{\dot{\rho }}}_{1,i}}}{{{\left( \delta _{2,i}^{2}\rho _{1,i}^{2}-{{s}^{2}} \right)}^{2}}}\text{d}s} \end{aligned} \right. $$ (52) 定义角速度跟踪误差为
$$ {{e}_{2,i}} = {{X}_{2,i}}-{{\chi }_{1,i}} $$ (53) 式中, $ {{\chi }_{1,i}} $为虚拟控制器, 将$ {{\chi }_{1,i}} $通过一阶滤波器, 可得
$$ {{\varepsilon }_{1,i}}{{\dot{\vartheta }}_{1,i}} = {{\chi }_{1,i}}-{{\vartheta }_{1,i}},\text{ }{{\vartheta }_{1,i}}(0) = {{\chi }_{1,i}}(0) $$ (54) 式中, $ {{\varepsilon }_{1,i}}\in {{\bf{R}}_{+}} $为滤波时间常数; $ {{\vartheta }_{1,i}} $和$ {{\dot{\vartheta }}_{1,i}} $分别为$ {{\chi }_{1,i}} $ 和 $ {{\dot{\chi }}_{1,i}} $ 的估计值, 定义其估计误差为 $ {{\zeta }_{1,i}} = {{\vartheta }_{1,i}}-{{\chi }_{1,i}} $, $ {{\xi }_{1,i}} = {{\dot{\vartheta }}_{1,i}}-{{\dot{\chi }}_{1,i}} $. 因此, $ {{\dot{V}}_{1}} $可改写为
$$ \begin{split} {{{\dot{V}}}_{1}} =\; &\sum\limits_{i = 1}^{3}{\left( {{\eta }_{1,i}}+{{\lambda }_{1,i}}({{X}_{2,i}}-{{{\dot{X}}}_{1d,i}}) \right)} = \\ & \sum\limits_{i = 1}^{3}{\left( {{\eta }_{1,i}}+{{\lambda }_{1,i}}({{e}_{2,i}}-{{{\dot{X}}}_{1d,i}}+{{\chi }_{1,i}}) \right)} = \\ & \sum\limits_{i = 1}^{3}{\left( {{\eta }_{1,i}}+{{\lambda }_{1,i}}({{e}_{2,i}}-{{{\dot{X}}}_{1d,i}}+{{\vartheta }_{1,i}}-{{\zeta }_{1,i}}) \right)} \end{split} $$ (55) 对式(55)应用Young不等式, 可得
$$ \begin{split} {{\lambda }_{1,i}}({{e}_{2,i}}-{{\zeta }_{1,i}})\le \;&\frac{2(1-{{\ell }_{1}})e_{1,i}^{2}\delta _{1,i}^{4}\rho _{1,i}^{4}}{{{\left( \delta _{1,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2} \right)}^{2}}}\;+ \\ {}& \frac{2{{\ell }_{1}}e_{1,i}^{2}\delta _{2,i}^{4}\rho _{1,i}^{4}}{{{\left( \delta _{2,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2} \right)}^{2}}}+\frac{e_{2,i}^{2}+\zeta _{1,i}^{2}}{2} = \\ {}&2\lambda _{1,i}^{2}+\frac{e_{2,i}^{2}+\zeta _{1,i}^{2}}{2} \\[-15pt]\end{split} $$ (56) 因此, 虚拟控制器$ {{\chi }_{1,i}} $可设计为
$$ \left\{ \begin{aligned} & {{\chi }_{1,i}} = {{{\dot{X}}}_{1d,i}}-2{{\lambda }_{1,i}}-\frac{{{\eta }_{1,i}}}{{{\lambda }_{1,i}}}-{{k}_{1,i}}{{e}_{1,i}},&{{e}_{1,i}}\ne 0 \\ & {{\chi }_{1,i}} = {{{\dot{X}}}_{1d,i}},&{{e}_{1,i}} = 0 \end{aligned} \right. $$ (57) 式中, $ {{k}_{1,i}}\in {{\bf{R}}_{+}} $为控制增益.
步骤 2. 考虑到$ {{\boldsymbol{e}}_{2}} = {{\boldsymbol{X}}_{2}}-{{\boldsymbol{\chi }}_{1}} $, 可得
$$ {{\dot{\boldsymbol{e}}}_{2}} = {{\dot{\boldsymbol{X}}}_{2}}-{{\dot{\boldsymbol{\chi }}}_{1}} = {{\boldsymbol{f}}_{1}}+{{{\boldsymbol{g}}}_{1}}{\boldsymbol{u}}+{{\boldsymbol{d}}_{1}}-{{\dot{\boldsymbol{\chi }}}_{1}} $$ (58) 为估计总扰动$ {{\boldsymbol{d}}_{1}} $, 可依据式(24)建立干扰观测器. 在步骤2中, 定义如下Lyapunov函数:
$$ {{V}_{2}} = {{V}_{1}}+\frac{1}{2}{\boldsymbol{e}}_{2}^{{\rm{T}}}{{\boldsymbol{e}}_{2}} $$ (59) 对式(59)求导可得
$$ {{\dot{V}}_{2}} = {{\dot{V}}_{1}}+{\boldsymbol{e}}_{2}^{{\rm{T}}}\left( {{\boldsymbol{f}}_{1}}+{{{\boldsymbol{g}}}_{1}}{\boldsymbol{u}}+{{\boldsymbol{d}}_{1}}-{{{\dot{\boldsymbol{\chi }}}}_{1}} \right) $$ (60) 考虑到$ {{\hat{\boldsymbol{d}}}_{1}} = {{\boldsymbol{z}}_{2}} $, $ {{{{\hat{\dot{\boldsymbol{\chi}} }}}}_{1}} = {{\dot{\vartheta }}_{1,i}} $, 姿态跟踪控制器可设计为
$$ \left\{ \begin{aligned} & {\boldsymbol{u}} = {\boldsymbol{g}}_{1}^{-1}{{{\boldsymbol{u}}}_{e}} \\ & {{u}_{e,i}} = -{{f}_{1,i}}-{{z}_{2,i}}+{{{\dot{\vartheta }}}_{1,i}}-\left({{k}_{2,i}}+\frac{3}{2}\right){{e}_{2,i}} \end{aligned} \right. $$ (61) 式中, $ {{k}_{2,i}}\in {{\bf{R}}_{+}} $为控制增益.
2.5 固定时间稳定性分析
本节对闭环系统进行稳定性分析, 各状态变量的有界性以及跟踪误差的收敛特性由定理2给出.
定理 2. 针对式(14)所描述的飞行器姿态跟踪控制系统, 基于固定时间干扰观测器(24)对总扰动$ {{\boldsymbol{d}}_{1}} $的精确估计, 设计控制器$ {\boldsymbol{u}} $(联立式(57)和式(61)), 可得到如下结论: 1) 闭环系统跟踪误差$ {{e}_{1,i}} $和$ {{e}_{2,i}} $均一致最终有界; 2) 姿态跟踪误差$ {{e}_{1,i}} $将在固定时间内收敛至预设性能边界(Prescribed performance bound, PPB)内; 3) 在整个时域内, 预先设定的性能约束条件(47)可得到满足.
证明. 联立式(55)和式(60), 可得
$$\begin{split} {{\dot{V}}_{2}} =\;& {{\dot{V}}_{1}}+\sum\limits_{i = 1}^{3}\bigg( -{{e}_{2,i}}{{e}_{z2,i}}+{{e}_{2,i}}{{\xi }_{1,i}}\;-\\ &\left({{k}_{2,i}}+\frac{3}{2}\right)e_{2,i}^{2} \bigg)\end{split} $$ (62) 对于式(62), 应用Young不等式, 可得
$$ \begin{split} & \sum\limits_{i = 1}^{3}{\left( -{{e}_{2,i}}{{e}_{z2,i}}+{{e}_{2,i}}{{\xi }_{1,i}} \right)}\le \\ &\qquad \sum\limits_{i = 1}^{3}{\left( \frac{e_{2,i}^{2}+e_{z2,i}^{2}}{2}+\frac{e_{2,i}^{2}+\xi _{1,i}^{2}}{2} \right)} \end{split} $$ (63) 对于$ {{\dot{V}}_{1}} $, 联立式(55) ~ (57), 存在如下关系:
$$ \begin{split} {{{\dot{V}}}_{1}}\le {}&-\sum\limits_{i = 1}^{3}\left[ {{k}_{1,i}}\rho _{1,i}^{2}e_{1,i}^{2}\left( \frac{(1-{{\ell }_{1}})\delta _{1,i}^{2}}{\delta _{1,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2}}\;+\right. \right. \\ {}& \left. \left. \frac{{{\ell }_{1}}\delta _{2,i}^{2}}{\delta _{2,i}^{2}\rho _{1,i}^{2}-e_{1,i}^{2}} \right) \right]+ \frac{1}{2}\sum\limits_{i = 1}^{3}{\left( e_{2,i}^{2}+\zeta _{1,i}^{2} \right)} \end{split} $$ (64) 将式(64)和式(63)代入式(62), 可得
$$ \begin{split} {{{\dot{V}}}_{2}}\le &-\sum\limits_{i = 1}^{3}{{{k}_{1,i}}{{\lambda }_{1,i}}{{e}_{1,i}}}+\frac{1}{2}\sum\limits_{i = 1}^{3}{\left( e_{2,i}^{2}+\zeta _{1,i}^{2} \right)}\;+ \\ & \sum\limits_{i = 1}^{3}{\left( e_{2,i}^{2}+\frac{e_{z2,i}^{2}+\xi _{1,i}^{2}}{2} \right)}+\\ &\sum\limits_{i = 1}^{3}{\left( -{{k}_{2,i}}e_{2,i}^{2}-\frac{3e_{2,i}^{2}}{2}\; \right)} = \\ &-\sum\limits_{i = 1}^{3}{{{k}_{1,i}}{{\lambda }_{1,i}}{{e}_{1,i}}}-\sum\limits_{i = 1}^{3}{{{k}_{2,i}}e_{2,i}^{2}}\;+ \\ & \frac{1}{2}\sum\limits_{i = 1}^{3}{\left( e_{z2,i}^{2}+\zeta _{1,i}^{2}+\xi _{1,i}^{2} \right)} \end{split} $$ (65) 根据文献[28]中定理1可知, 对于式(49)所示积分障碍Lyapunov函数, 当$ {{e}_{1,i}} $始终处于式(47)所述性能约束之内时, $ {{V}_{B,i}}({{e}_{1,i}},{{\rho }_{1,i}})\le {{\lambda }_{1,i}}{{e}_{1,i}} $恒成立. 因此, 可得到如下关系:
$$ \begin{split} & -\sum\limits_{i = 1}^{3}{{{k}_{1,i}}{{\lambda }_{1,i}}{{e}_{1,i}}}-\sum\limits_{i = 1}^{3}{{{k}_{2,i}}e_{2,i}^{2}} \le\\ &\qquad -\sum\limits_{i = 1}^{3}{{{k}_{1,i}}{{V}_{B,i}}({{e}_{1,i}},{{\rho }_{1,i}})}-\sum\limits_{i = 1}^{3}{{{k}_{2,i}}e_{2,i}^{2}} \le\\ &\qquad -k_{1}^{*}\sum\limits_{i = 1}^{3}{{{V}_{B,i}}({{e}_{1,i}},{{\rho }_{1,i}})}-2k_{2}^{*}\sum\limits_{i = 1}^{3}{{e_{2,i}^{2}}/{2}\;} = \\ &\qquad -k_{1}^{*}{{V}_{1}}-2k_{2}^{*}({{V}_{2}}-{{V}_{1}}) = (2k_{2}^{*}-k_{1}^{*}){{V}_{1}}-2k_{2}^{*}{{V}_{2}} \end{split} $$ (66) 式中, $ k_{1}^{*} = \min \left\{ {{k}_{1,i}} \right\} $, $ k_{2}^{*} = \min \left\{ {{k}_{2,i}} \right\} $, $ i = 1,2,3 $. 然后, 基于文献[29]和第2.2节所述固定时间干扰观测器收敛特性, 干扰观测器观测误差$ {{e}_{z2,i}} $, 滤波器估计误差$ {{\zeta }_{1,i}} $, $ {{\xi }_{1,i}} $均有界, 并且均可通过参数整定收敛至极小值. 换言之, 总存在如下关系:
$$ \frac{1}{2}\sum\limits_{i = 1}^{3}{\left( e_{z2,i}^{2}+\zeta _{1,i}^{2}+\xi _{1,i}^{2} \right)}\le {{\imath }_{1}} $$ (67) 式中, $ {{\imath}_{1}}\in {{\bf{R}}_{+}} $为一足够小的常数. 因此, $ {{\dot{V}}_{2}} $满足:
$$ \begin{split} {{{\dot{V}}}_{2}}\le {}& (2k_{2}^{*}-k_{1}^{*}){{V}_{1}}-2k_{2}^{*}{{V}_{2}}\;+ \\ {}& \frac{1}{2}\sum\limits_{i = 1}^{3}{\left( e_{z2,i}^{2}+\zeta _{1,i}^{2}+\xi _{1,i}^{2} \right)} \le -{{\hbar }_{1}}{{V}_{2}}+{{\imath}_{1}} \end{split} $$ (68) 式中, $ {{\hbar }_{1}} = \min \left\{ 2k_{2}^{*},k_{1}^{*} \right\} $. 解不等式(68), 可得
$$ \begin{split} 0\le \;&{{V}_{2}}(t)\le \left( V(0)-\frac{{{\imath}_{1}}}{{{\hbar }_{1}}} \right){{e}^{-{{\hbar }_{1}}t}}+\frac{{{\imath}_{1}}}{{{\hbar }_{1}}} \le\\ {}& {{V}_{2}}(0){{e}^{-{{\hbar }_{1}}t}}+\frac{{{\imath}_{1}}}{{{\hbar }_{1}}} \end{split} $$ (69) 基于式(69)和引理1, 可知闭环系统跟踪误差$ {{e}_{1,i}} $和$ {{e}_{2,i}} $均一致最终有界. 进一步根据定义2和引理2可知, 姿态跟踪误差$ {{e}_{1,i}} $始终不超出预设性能约束(47). 与此同时, 当$ {{e}_{1,i}} $始终满足性能约束时, 基于式(47)所建立的不等式关系, 易证明$ {{e}_{1,i}} $将在固定时间内收敛至预设的极小边界内.
□ 注 6. 本文闭环控制系统的固定时间稳定性由所设计性能函数保证, 且该性能函数具有如下特性: 系统调节时间上界与控制参数和初始状态无关, 而是根据实际需要预先设定. 根据文献[30]可知, 当稳态边界值$ {{\rho }_{\infty }} $足够小时, 可认为跟踪误差系统是固定时间稳定的.
注 7. 与现有研究中的有限时间或固定时间控制方法(如文献[2]和[4])相比, 本节方法避免了状态量的分数次幂和符号函数的使用, 极大地降低了计算量和简化了控制器结构. 此外, 固定时间干扰观测器可由线性观测器代替, 仍可基于性能函数获得优异的快速收敛能力, 实际中可以根据计算能力和收敛性能的实际需求选择相应的观测器. 这也进一步说明了所提出的控制方法的灵活性.
3. 仿真验证与结果分析
本节将基于所建立飞行器运动模型和所设计姿态控制器开展仿真验证. 针对飞行器滑翔飞行过程中的某一段, 进行姿态控制基本性能仿真, 以验证飞行器姿态跟踪性能和所提出控制方法的有效性; 同时, 为了进一步验证所提出方法的鲁棒性, 进行参数拉偏和注入干扰等相关内容的适应性仿真.
仿真1: 基本性能仿真. 在不考虑模型不确定性和外部干扰的情况下, 选取某一工作点进行飞行器滑翔段六自由度运动仿真, 飞行器机体参数如表2所示, 仿真初始条件设置如表3所示. 具体来说:
表 2 高超声速变外形飞行器机体参数Table 2 Body parameters of HMV参量 符号 数值 单位 机身质量 $m_f$ 2950 kg 折叠翼质量 $m_1, m_2$ 55 kg $x$主轴转动惯量 $I_{xx}$ $\left [ 283,298 \right ] $ kg·m2 $y$主轴转动惯量 $I_{yy}$ $\left [ 2\;679,2\;722 \right ]$ kg·m2 $z$主轴转动惯量 $I_{zz}$ $\left [ 2\;528,2\;630 \right ]$ kg·m2 惯量积 $I_{xy}$ $\left [ 163,169 \right ] $ kg·m2 参考面积 $S_r$ 1.8 m2 参考气动弦长 $c_A$ 2.4 m 参考气动展长 $b_A$ 1.1 m 表 3 仿真参数设置Table 3 Setting of simulation parameters参数类型 参数值 初始状态参数 $H=35$ km, $V=3\;200$ m/s $\lambda ={{120}^{\circ }}$, $\phi ={{20}^{\circ}}$, $\theta=-{{1}^{\circ}}$, ${{\psi}_{v}}={{10}^{\circ}}$ $\alpha={{8}^{\circ}}$, $\beta={{1}^{\circ}}$, $\sigma={{18}^{\circ}}$ ${{\omega}_{x}}={{\omega}_{y}}={{\omega}_{z}}=0$, ${{\delta}_{x}}={{\delta}_{y}}={{\delta}_{z}}=0$ 控制参数 ${{\boldsymbol{\rho }}_{0}}={{\left[ {{\rho }_{0,1}},{{\rho }_{0,2}},{{\rho }_{0,3}} \right]}^{\text{T}}}={{\left[ 5,3,5 \right]}^{\text{T}}}$ ${{\boldsymbol{\rho }}_{\infty }}={{\left[ {{\rho }_{\infty ,1}},{{\rho }_{\infty ,2}},{{\rho }_{\infty ,3}} \right]}^{\text{T}}}={{\left[ 0.2,0.1,0.3 \right]}^{\text{T}}}$ ${{m}_{1,i}}=3$, ${{m}_{2,i}}=5$, ${{n}_{1,i}}=5$, ${{n}_{2,i}}=7$ ${{\alpha }_{01,i}}=0.15$, ${{\alpha }_{02,i}}=0.2$, ${{\delta }_{1,i}}={{\delta }_{2,i}}=1$ ${{k}_{1,i}}={{k}_{2,i}}=2$, $\text{ }{{k}_{3,i}}={{k}_{4,i}}=4$ ${{\varepsilon }_{1,i}}=0.02$, ${{\gamma }_{1,i}}=0.6$, ${{\gamma }_{2,i}}=1.4$ ${{k}_{z1,i}}=4$, ${{k}_{z2,i}}=4$, ${{\varepsilon }_{0,i}}=0.2$ 仿真步长 d$t$= 0.01 s 外部干扰项 $\Delta {{d}_{1,1}}=500\left( -\cos ({\pi t}/{20})+\sin ({\pi t}/{40}) \right)\;\text{N}\cdot \text{m}$ $\Delta {{d}_{1,2}}=300\left( -\cos ({\pi t}/{30})+\sin ({\pi t}/{60}) \right)\;\text{N}\cdot \text{m}$ $\Delta {{d}_{1,3}}=1\;000\cos ({\pi t}/{30})\sin ({\pi t}/{20})\;\text{N}\cdot \text{m}$ 模型不确定项 $\Delta{{C}_{L}}=\Delta{{C}_{D}}=\Delta{{C}_{Y}}=\pm20\%$ $\Delta{{C}_{mx}}=\Delta{{C}_{my}}=\Delta{{C}_{mz}}=\pm20\%$ $\Delta{{I}_{xx}}=\Delta{{I}_{yy}}=\Delta{{I}_{zz}}=\Delta{{I}_{xy}}=\pm20\%$ $\Delta{{S}_{r}}=\Delta{{b}_{A}}=\Delta{{c}_{A}}=\pm5\%$ 1) 采用式(57)和式(61)所示控制器, 不采用观测器, 参数设置如表3所示.
2) 为验证飞行器在不同变形条件(变形量和变形速率)下的姿态跟踪性能和动态特性, 分别设计了多段不同的折叠角剖面.
3) 姿态角期望指令为预先设定, 如图10所示.
仿真2: 适应性仿真. 考虑模型不确定性和外部干扰, 总扰动参数设置如表3所示. 选取与仿真1相同的初始条件, 进行飞行器滑翔段六自由度运动仿真. 具体来说:
1) 采用式(57)和式(61)所示控制器, 并采用式(24)所示干扰观测器对总扰动进行精确估计, 参数设置如表3所示.
2) 飞行器变外形策略一般由前端轨迹及制导回路给出, 在仿真中, 折叠角和期望姿态角均为预先给定的开环剖面.
3) 将本文所提出控制方法(记为FxTPPC) 分别与文献[31]中的动态面控制方法(记为DSC)和文献[32]中的预设性能动态面控制方法(记为PPDSC)进行对比, 且均采用式(24)所示干扰观测器.
在仿真中, 为定量描述3种方法下的姿态跟踪性能, 定义姿态跟踪累计误差CE, CE = $ \int_{0}^{T_f} t{\left \| {\boldsymbol{e}}_1 \right \| }\text{d}t $. 对于第1.4节所述外部干扰和模型不确定性, 本文仿真通过干扰注入和参数拉偏进行模拟, 具体见表3.
图10 ~ 12分别给出了仿真1的姿态角跟踪曲线、姿态角跟踪误差曲线和姿态角速度变化曲线. 由图10 ~ 12可知, 可将整个仿真过程分为3段: 1) 0 ~ 4 s, 姿态角跟踪误差从一较大值收敛至零, 该过程中飞行器姿态角速度剧烈变化; 2) 4 ~ 10 s, 姿态角跟踪误差已收敛至零附近, 姿态角速度也趋近于零; 3) 10 ~ 20 s, 跟踪一组三角函数型姿态角指令, 姿态角速度持续变化.
由图13可知, 在每一段, 均令折叠角分别以$ 2^{\circ}/{\rm{s}} $, $ 6^{\circ}/{\rm{s}}$和$ 10^{\circ} /{\rm{s}}$的速度变化, 直至区间终点. 结合图10、11可知, 在不同变形条件, 采用本文所提出的控制方法均可以实现较高品质的控制性能. 在0 ~ 4 s阶段, 各通道稳态误差均能在预设时间2 s内收敛至$ {10^{-3}}^{\circ} $内, 并在4 ~ 10 s阶段保持这一误差区间; 在10 ~ 20 s阶段, 攻角误差不超过$ 0.05^{\circ} $, 侧滑角误差不超过$ 2\times {10^{-4}}^{\circ} $, 倾侧角误差不超过$ 0.1^{\circ} $.
图14 ~ 19分别给出了仿真中各段的气动力/力矩和附加力/力矩变化曲线. 综合图14 ~ 19, 可得到如下结论:
图 14 仿真1气动力和附加力变化曲线 (0 ~ 3 s)Fig. 14 Curves of aerodynamic force and additional force in Simulation 1 (0 ~ 3 s)
图 15 仿真1气动力和附加力变化曲线 (5 ~ 9 s)Fig. 15 Curves of aerodynamic force and additional force in Simulation 1 (5 ~ 9 s)
图 16 仿真1气动力和附加力变化曲线 (11 ~ 16 s)Fig. 16 Curves of aerodynamic force and additional force in Simulation 1 (11 ~ 16 s)
图 17 仿真1气动力矩和附加力矩变化曲线 (0 ~ 3 s)Fig. 17 Curves of aerodynamic torque and additional torque in Simulation 1 (0 ~ 3 s)
图 18 仿真1气动力矩和附加力矩变化曲线 (5 ~ 9 s)Fig. 18 Curves of aerodynamic torque and additional torque in Simulation 1 (5 ~ 9 s)
图 19 仿真1气动力矩和附加力矩变化曲线 (11 ~ 16 s)Fig. 19 Curves of aerodynamic torque and additional torque in Simulation 1 (11 ~ 16 s)1) 机翼折叠变形对气动力影响较大, 主要是对飞行器升阻特性的影响, 但变形对气动力矩影响较小. 机翼折叠变形对附加力矩影响远大于对附加力的影响. 总体来看, 折叠角变化率越大, 所产生的附加力及附加力矩越大.
2) 折叠变形引起的附加力与气动力相差较大, 附加力远小于相应通道的阻力、升力和侧力. 折叠变形引起的附加力矩与气动力矩基本相当, 并且, 俯仰通道的附加力矩最大, 远高于偏航通道附加力矩, 滚转通道附加力矩最小. 由此说明变形机构产生的附加力矩是飞行器所受力矩的主要来源, 执行机构需产生相应控制力矩抵消附加力矩并达到跟踪期望姿态指令的目标, 也进一步说明了本文控制方法对变外形的适应能力.
3) 变形过程中, 附加力和附加力矩受姿态角速度变化影响较大, 远大于其受折叠角变化率的影响. 尤其是初始姿态稳定控制阶段, 当姿态角速度较大时, 变形所导致的附加力及附加力矩将急剧增大. 这说明附加力及附加力矩与状态(尤其是姿态角速度)严重耦合, 若要实现快速收敛的瞬态性能, 必然引起极大的附加力矩, 但本文控制方法仍可以实现较高品质的控制性能.
图20 ~ 22分别给出了仿真2攻角、侧滑角和倾侧角的跟踪曲线; 图23给出了仿真2总扰动及其观测误差的变化曲线; 图24 ~ 26分别给出了仿真2等效舵偏角、折叠角和姿态角跟踪累积误差的变化曲线. 结合图20 ~ 22和图26可知, 在收敛时间上, FxTPPC方法最短, PPDSC方法次之, DSC方法最长; 在稳态精度上, FxTPPC方法最优, 可确保姿态角跟踪误差始终处于预设性能约束内, 而DSC方法则不能全程保证, 且稳态误差最大. 从图23可以看出, 本文所设计的固定时间干扰观测器可实现对模型不确定性和外部干扰共同构成的复杂总扰动的精确估计. 由图24可知, 在FxTPPC方法下, 舵偏角平滑变化, 无饱和, 无抖振. 综合来看, 即使在模型不确定性和外部干扰作用下, 本文提出的FxTPPC方法仍能实现较高品质的姿态控制, 相比于PPDSC方法和DSC方法, FxTPPC方法具备更优异的控制性能, 即收敛速度更快, 超调量更小, 稳态误差更小, 进一步说明了本文方法的有效性与鲁棒性.
图 23 仿真2总扰动及其观测误差曲线Fig. 23 Curves of total disturbance and its observation error in Simulation 24. 结论
针对一种新型折叠式高超声速变外形飞行器, 本文建立了飞行器完整运动模型并开展了气动特性分析和运动特性分析, 在此基础上, 充分考虑不同变形条件、模型不确定性和外部干扰对飞行器姿态控制的影响, 融合基于固定时间干扰观测器的扰动补偿控制、基于固定时间性能函数的预设性能控制和动态面控制的复合控制方案, 设计了一种新型固定时间姿态控制方法. 此外, 本文还进行了滑翔段某工况下的对比仿真, 仿真结果表明, 本文所提出的控制方法在大尺度变形、模型不确定性和外部干扰影响下可实现预先设定的姿态跟踪瞬态和稳态性能, 充分说明了本文所提出方法的有效性和鲁棒性. 在今后的工作中, 为进一步提高高超声速变外形飞行器在高不确定性和强干扰环境下的姿态控制精度和控制器的鲁棒性, 可从以下几个方面进一步完善本文工作: 1) 如何降低初值、过程噪声等对观测器的影响, 并进一步提升观测器系统的收敛能力, 是进一步提升本文控制方法鲁棒性的关键; 2) 实际中往往存在执行机构控制能力不足的问题, 如何在保证姿态跟踪瞬态性能和稳态性能的同时降低本文控制方法在姿态响应初始时刻对控制量的需求以避免控制饱和, 也是后续需要关注的重点.
-
图 1 商业住宅区域综合能源系统架构
Fig. 1 Integrated energy system architecture for commercial and residential area
图 5 日均冷负荷与光伏预测功率曲线
Fig. 5 Average daily cooling load and photovoltaic predicted power curves
图 6 日均电负荷与日均热负荷曲线
Fig. 6 Daily average electric load and daily average heat load curve
图 9 优化前后内燃机出力对比
Fig. 9 Comparison of internal combustion engines before and after optimize output
表 1 收敛度对比
Table 1 Convergence contrast
算法 指标 ZDT1 ZDT2 ZDT3 AMOWOA M 9.41${\times{10^{-4}}}$ 9.59${\times{10^{-4}}}$ 9.68${\times{10^{-4}}}$ V 2.26${\times{10^{-5}}}$ 3.41${\times{10^{-5}}}$ 2.16${\times{10^{-5}}}$ NSGA-II M 9.79${\times{10^{-4}}}$ 9.68${\times{10^{-4}}}$ 9.84${\times{10^{-4}}}$ V 4.88${\times{10^{-5}}}$ 5.84${\times{10^{-5}}}$ 3.63${\times{10^{-5}}}$ MOPSO M 9.46${\times{10^{-4}}}$ 1.42${\times{10^{-3}}}$ 9.73${\times{10^{-4}}}$ V 3.42${\times{10^{-5}}}$ 8.26${\times{10^{-5}}}$ 3.79${\times{10^{-5}}}$ PESA-II M 1.05${\times{10^{-3}}}$ 7.40${\times{10^{-4}}}$ 7.89${\times{10^{-3}}}$ V 0.00 0.00 1.10${\times{10^{-4}}}$ NSPSO M 6.42${\times{10^{-3}}}$ 9.51${\times{10^{-3}}}$ 4.91${\times{10^{-3}}}$ V 0.00 0.00 0.00 
下载: 导出CSV
表 2 多样度对比
Table 2 Diversity contrast
算法 指标 ZDT1 ZDT2 ZDT3 AMOWOA M 0.65560 0.74680 0.79080 V 0.02109 0.03116 0.02679 NSGA-II M 0.74470 0.87290 0.78760 V 0.02901 0.05793 0.06771 MOPSO M 0.75250 0.93860 0.95170 V 0.03574 0.06475 0.01563 PESA-II M 0.84810 0.89290 1.22730 V 0.00287 0.05740 0.02930 NSPSO M 0.90700 0.92200 0.06210 V 0.00 1.20${\times{10^{-4}}}$ 6.90${\times{10^{-4}}}$
下载: 导出CSV
表 3 设备规格
Table 3 Specification of equipment
设备 配置容量 能效系数 (COP) 内燃机 10 000 kW — 光伏 7 100 kW — 电制冷机 2 000 kW 3.1 热泵 5 000 kW 4.4 (热)/5 (冷) 溴化锂余 8 000 kW 1.0 热机组 蓄电池 6 000 kWh 0.9 (充/放) 储热设备 5 000 kWh 0.9 (充/放) 储冷设备 2 000 kWh 0.9 (充/放)
下载: 导出CSV
表 4 模型参数
Table 4 Model parameter
参数 数值 内燃机电效率 41.33% 内燃机热效率 40.54% 内燃机燃料热耗率 7 962.726 kJ/kWh 电网输电效率 92% 发电厂发电效率 37%
下载: 导出CSV
表 5 初始运行条件
Table 5 Initial operating conditions
时段 内燃机出力 (kW) 1 (0:00−4:00) 4 000 2 (4:00−8:00) 4 000 3 (8:00−12:00) 8 000 4 (12:00−16:00) 8 000 5 (16:00−20:00) 8 000 6 (20:00−24:00) 4 000
下载: 导出CSV
A1 多目标优化标准测试函数表达式
A1 Multi-objective optimization standard test functions expression
测试函数 表达式 ZDT1 $\left\{\begin{aligned} &\min{f}_{1}\left({x}_{1}\right)={x}_{1}\\& \mathrm{min}{f}_{2}\left(x\right)=g\left(1-\sqrt{\frac{ {f}_{1} }{g}}\right)\\ &g\left(x\right)=1 +\frac{9\sum\limits _{i=2}^{m}{x}_{i} }{m-1}\end{aligned}\right.$ ZDT2 $\left\{\begin{aligned} &\min{f}_{1}\left({x}_{1}\right)={x}_{1}\\& \mathrm{min}{f}_{2}\left(x\right)=g\left(1-{\left(\frac{ {f}_{1} }{g}\right)}^{2}\right)\\& g\left(x\right)=1 +\frac{9\sum\limits _{i=2}^{m}{x}_{i} }{m-1}\end{aligned}\right.$ ZDT3 $\left\{\begin{aligned} &\min{f}_{1}\left({x}_{1}\right)={x}_{1}\\& \mathrm{min}{f}_{2}\left(x\right)=g\left(1-\sqrt{ \frac{ {f}_{1} }{g} }-\left(\frac{ {f}_{1} }{g}\right)\mathrm{sin}\left(10\pi {f}_{1}\right)\right)\\& g\left(x\right)=1 +\frac{9\sum\limits _{i=2}^{m}{x}_{i} }{m-1}\end{aligned}\right.$
下载: 导出CSV
-
[1] 王伟亮, 王丹, 贾宏杰, 陈沼宇, 郭炳庆, 周海明, 范孟华.能源互联网背景下的典型区域综合能源系统稳态分析研究综述[J]. 中国电机工程学报, 2016, 36(12):3292-3306Wang Wei-Liang, Wang Dan, Jia Hong-Jie, Chen Zhao-Yu, Guo Bing-Qing, Zhou Hai-Ming, Fan Meng-Hua.Review on steady state analysis of typical regional integrated energy system under the background of energy Internet[J]. Chinese Journal of Electrical Engineering, 2016, 36 (12): 3292-3306 [2] Kim I, James J A, Crittenden J. The case study of combined cooling heat and power and photovoltaic systems for building customers using HOMER software[J]. Electric Power Systems Research, 2017, 143:490-502 doi: 10.1016/j.jpgr.2016.10.061 [3] Gu W, Wang Z H, Wu Z, Luo Z, Tang Y Y, Wang J. An online optimal dispatch schedule for CCHP microgrids based on model predictive control[J]. IEEE Transactions on Smart Grid, 2017, 8 (5):2332-2342 [4] Zhang X, Shahidehpour M, Alabdulwahab A, Abusorrah A. Optimal expansion planning of energy hub with multiple energy infrastructures[J]. IEEE Transactions on Smart Grid, 2015, 6 (5):2302-2311 doi: 10.1109/TSG.2015.2390640 [5] 瞿凯平, 黄琳妮, 余涛, 张孝顺.碳交易机制下多区域综合能源系统的分散调度[J]. 中国电机工程学报, 2018, 38(03):697-707Qu Kai-Ping, Huang Lin-Ni, Yu Tao, Zhang Xiao-Shun. Decentralized scheduling of multi regional integrated energy system under carbon trading mechanism[J]. Chinese Journal of Electrical Engineering, 2018, 38 (03): 697-707 [6] 吴克河, 王继业, 李为, 朱亚运.面向能源互联网的新一代电力系统运行模式研究[J]. 中国电机工程学报, 2019, 39(04):966-979Wu Ke-He, Wang Ji-Ye, Li Wei, Zhu Ya-Yun. Research on operation mode of new generation power system for energy Internet[J]. Chinese Journal of Electrical Engineering, 2019, 39 (04): 966-979 [7] 别朝红, 林超凡, 李更丰, 邱爱慈.能源转型下弹性电力系统的发展与展望[J]. 中国电机工程学报, 2020, 40(09):2735-2745Bie Chao-Hong, Lin Chao-Fan, Li Geng-Feng, Qiu Ai-Ci. Development and prospect of flexible power system under energy transformation[J]. Chinese Journal of Electrical Engineering, 2020, 40 (09): 2735-2745 [8] Wang Y R, Zeng B, Guo J, Shi J Q, Zhang J H. Multi energy flow calculation method of electric heat gas integrated energy system[J]. Power Grid Technology, 2016, 40 (10): 2942-2951 [9] Hajabdollahi H, Ganjehkaviri A, Jaafar M N M. Assessment of new operational strategy in optimization of CCHP plant for different climates using evolutionary algorithms[J]. Applied Thermal Engineering, 2015, 75:468-480 doi: 10.1016/j.applthermaleng.2014.09.033 [10] 刘涤尘, 马恒瑞, 王波, 高文忠, 王骏, 闫秉科. 含冷热电联供及储能的区域综合能源系统运行优化. 电力系统自动化, 2018, 42(4): 113−120, 141Liu Di-Chen, Ma Heng-Rui, Wang Bo, Gao Wen-Zhong, Wang Jun, Yan Bing-Ke. Operation optimization of regional integrated energy system with combined cooling, heating and power supply and energy storage. Power System Automation, 2018, 42(4): 113−120, 141 [11] 施泉生, 丁建勇, 刘坤, 晏伟. 含电、气、热 3 种储能的微网综合能源系统经济优化运行. 电力自动化设备, 2019, 39(8): 269−276, 293Shi Quan-Sheng, Ding Jian-Yong, Liu Kun, Yan Wei. Economic optimization of microgrid integrated energy system with three kinds of energy storage including electricity, gas and heat. Power Automation Equipment, 2019, 39(8): 269−276, 293 [12] 顾洁, 白凯峰, 时亚军.基于多主体主从博弈优化交互机制的区域综合能源系统优化运行[J]. 电网技术, 2019, 43(09):3119-3134Gu Jie, Bai Kai-Feng, Shi Ya-Jun. Optimal operation of regional comprehensive energy system based on multi-agent master-slave game optimization interaction mechanism[J]. Power Grid Technology, 2019, 43 (09): 3119-3134 [13] 郑亚锋, 魏振华, 王春雨.计及储热装置的综合能源系统分层优化调度[J]. 中国电机工程学报, 2019, 39(S1):36-43 doi: 10.13334/J.0258-8013.PCSEE.190221Zheng Ya-Feng, Wei Zhen-Hua, Wang Chun-Yu. Hierarchical optimal scheduling of integrated energy system considering heat storage device[J]. Chinese Journal of Electrical Engineering, 2019, 39 (S1): 36-43 doi: 10.13334/J.0258-8013.PCSEE.190221 [14] 方彤, 蒋东方, 杨洋, 袁铁江.基于NSGA-II和熵权法的氢综合能源系统商业运营模式[J]. 中国电力, 2022, 55(01):110-118Fang Tong, Jiang Dong-Fang, Yang Yang, Yuan Tie-Jiang. Commercial operation mode of hydrogen integrated energy system based on NSGA-II and entropy weight method[J]. China Power, 2022, 55 (01): 110-118 [15] 齐世雄, 王秀丽, 邵成成, 王智冬, 朱承治.计及弹性恢复的区域综合能源系统多目标优化调度[J]. 中国电力, 2019, 52(06):19-26Qi Shi-Xiong, Wang Xiu-Li, Shao Cheng-Cheng, Wang Zhi-Dong, Zhu Cheng-Zhi.Multi objective optimal scheduling of regional integrated energy system considering elastic recovery[J]. China Power, 2019, 52 (06): 19-26 [16] 王磊, 姜涛, 宋丹, 崔杨, 胡扬宇, 柴旭峥, 唐耀华.基于灵活热电比的区域综合能源系统多目标优化调度[J]. 电力系统保护与控制, 2021, 49(08):151-159 doi: 10.19783/j.cnki.pspc.201561Wang Lei, Jiang Tao, Song Dan, Cui Yang, Hu Yang-Yu, Chai Xu-Zheng, Tang Yao-Hua. Multi objective optimal scheduling of regional integrated energy system based on flexible thermal power ratio[J]. Power System Protection and Control, 2021, 49 (08): 151-159 doi: 10.19783/j.cnki.pspc.201561 [17] 华煌圣, 刘育权, 熊文, 徐航, 施云辉, 董树锋.考虑综合能效水平的能源系统多目标优化运行[J]. 南方电网技术, 2018, 12(03):81-84 doi: 10.13648/j.cnki.issn1674-0629.2018.03.011Hua Huang-Sheng, Liu Yu-Quan, Xiong Wen, Xu Hang, Shi Yun-Hui, Dong Shu-Feng. Multi objective optimal operation of energy system considering comprehensive energy efficiency level[J]. China Southern Power Grid Technology, 2018, 12 (03): 81-84 doi: 10.13648/j.cnki.issn1674-0629.2018.03.011 [18] 程亮.基于鲁棒优化方法的区域综合能源系统运行优化调度策略[J]. 电力系统装备, 2021(18):2Cheng Liang.Operation optimization scheduling strategy of regional integrated energy system based on robust optimization method[J]. Power System Equipment, 2021 (18): 2 [19] 张涛, 郭玥彤, 李逸鸿, 余利, 章佳莹. 计及电气热综合需求响应的区域综合能源系统优化调度[J]. 电力系统保护与控制, 2021, 49(1):10 doi: 10.19783/j.cnki.pspc.200167Zhang Tao, Guo Yue-Tong, Li Yi-Hong, Yu Li, Zhang Jia-Ying. Optimal dispatching of regional integrated energy system considering comprehensive demand response of electrical and thermal[J]. Power System Protection and Control, 2021, 49 (1): 10 doi: 10.19783/j.cnki.pspc.200167 [20] 魏震波, 任小林, 黄宇涵. 考虑综合需求侧响应的区域综合能源系统多目标优化调度[J]. 电力建设, 2020, 41(7):8 doi: 10.12204/j.issn.1000-7229.2020.07.012Wei Zhen-Bo, Ren Xiao-Lin, Huang Yu-Han. Multi objective optimal scheduling of regional integrated energy system considering comprehensive demand side response[J]. Power Construction, 2020, 41 (7): 8 doi: 10.12204/j.issn.1000-7229.2020.07.012 [21] 耿琪, 胡炎, 何建宗, 周永言, 赵伟. 基于纳什谈判的区域综合能源系统运行优化[J]. 电力建设, 2020, 41(1):114-125Geng Qi, Hu Yan, He Jian-Zong, Zhou Yong-Yan, Zhao Wei.Operation optimization of regional integrated energy system based on nash negotiation[J]. Power Construction, 2020, 41 (1): 114-125 [22] 施云辉, 郭创新. 考虑运行风险的含储能综合能源系统优化调度[J]. 发电技术, 2020, 41(1):8 doi: 10.12096/j.2096-4528.pgt.19142Shi Yunhui, Guo Chuang-Xin. Optimal scheduling of integrated energy system with energy storage considering operation risk[J]. Power Generation Technology, 2020, 41 (1): 8 doi: 10.12096/j.2096-4528.pgt.19142 [23] 何畅, 程杉, 徐建宇, 等. 基于多时间尺度和多源储能的综合能源系统能量协调优化调度[J]. 电力系统及其自动化学报, 2020, 32(2):9 doi: 10.19635/j.cnki.csu-epsa.000363He Chang, Cheng Shan, Xu Jian-Yu, et al. Energy coordinated optimal scheduling of integrated energy system based on multi time scale and multi-source energy storage[J]. Journal of Power System and Automation, 2020, 32 (2): 9 doi: 10.19635/j.cnki.csu-epsa.000363 [24] 崔杨、闫石、仲悟之、王铮、张鹏、赵钰婷. 含电转气的区域综合能源系统热电优化调度[J]. 电网技术, 2020, 44(11):10 doi: 10.13335/j.1000-3673.pst.2019.2468Cui Yang, Yan Shi, Zhong Wu-Zhi, Wang Zheng, Zhang Peng, Zhao Yu-Ting. Thermoelectric optimal scheduling of regional integrated energy system with electricity to gas[J]. Power Grid Technology, 2020, 44 (11): 10 doi: 10.13335/j.1000-3673.pst.2019.2468 [25] 杨雍琦, 薛万磊, 张海静, 等. 计及需求响应的区域综合能源系统的优化调度方法[J]. 中国电力, 2021, 54(4):10Yang Yong-Qi, Xue Wan-Lei, Zhang Hai-Jing, et al.Optimal scheduling method of regional integrated energy system considering demand response[J]. China Power, 2021, 54 (4): 10 [26] 李玉帅, 李天义, 高炜, 高文忠. 基于异步动态事件触发通信策略的综合能源系统分布式协同优化运行方法. 自动化学报, 2020, 46(9): 1831-1843 doi: 10.16383/j.aas.c200172Li Yu-Shuai, Li Tian-Yi, Gao Wei, Gao Wen-Zhong. Distributed cooperative optimal operation method of integrated energy system based on asynchronous dynamic event triggered communication strategy. Journal of Automation, 2020, 46 (9): 1831 - 1843 doi: 10.16383/j.aas.c200172 [27] 耿志强, 毕帅, 王尊, 朱群雄, 韩永明. 基于改进NSGA-II算法的乙烯裂解炉操作优化[J]. 化工学报, 2020, 71(3):7Geng Zhi-Qiang, Bi Shuai, Wang Zun, Zhu Qun-Xiong, Han Yong-Ming. Operation optimization of ethylene cracking furnace based on improved NSGA - II algorithm[J]. Journal of Chemical Engineering, 2020, 71 (3): 7 [28] 曾鸣, 韩旭, 李源非, 刘金洁, 彭丽霖.基于Tent映射混沌优化NSGA-II算法的综合能源系统多目标协同优化运行[J]. 电力自动化设备, 2017, 37(06):220-228Zeng Ming, Han Xu, Li Yuan-Fei, Liu Jin-Jie, Peng Li-Lin. Multi objective cooperative optimization operation of integrated energy system based on tent map chaotic optimization NSGA - II algorithm[J]. Power Automation Equipment, 2017, 37 (06): 220-228 [29] 董帅, 王成福, 梁军, 董晓明, 梁正堂.等.计及电转气运行成本的综合能源系统多目标日前优化调度[J]. 电力系统自动化, 2018, 42(11):8-15 doi: 10.7500/AEPS20170721003Dong Shuai, Wang Cheng-Fu, Liang Jun, Dong Xiao-Ming, Liang Zheng-Tang.Multi objective day ahead optimal scheduling of integrated energy system considering the operation cost of electricity to gas[J]. Power System Automation, 2018, 42 (11): 8-15 doi: 10.7500/AEPS20170721003 [30] Mirjalili S, Lewis A. The whale optimization algorithm[J]. Advances in Engineering Software, 2016, 95: 51-67 doi: 10.1016/j.advengsoft.2016.01.008 [31] Kaur G, Arora S. Chaotic whale optimization algorithm[J]. Journal of Computational Design and Engineering, 2018, 5:275–284 doi: 10.1016/j.jcde.2017.12.006 [32] 黄博南, 王勇, 李玉帅, 刘鑫蕊, 杨超. 基于分布式神经动态优化的综合能源系统多目标优化调度[J]. 自动化学报, 2020, 46: 1-19 doi: 10.16383/j.aas.c200168Huang Bo-Nan, Wang Yong, Li Yu-Shuai, Liu Xin-Rui, Yang Chao.Multi objective optimal scheduling of integrated energy system based on distributed neural dynamic optimization[J]. Journal of Automation, 2020, 46: 1–19 doi: 10.16383/j.aas.c200168 [33] 孙才志, 林学钰. 基于层次分析的模糊一致性判断矩阵及其应用[J]. 模糊系统与数学, 2002, 16(3): 59-63 doi: 10.3969/j.issn.1001-7402.2002.03.011Sun Cai-Zhi, Lin Xue-Yu.Fuzzy consistency judgment matrix based on analytic hierarchy process and its application[J]. Fuzzy Systems and Mathematics, 2002, 16 (3): 59-63 doi: 10.3969/j.issn.1001-7402.2002.03.011 [34] 陈民铀, 张聪誉, 罗辞勇. 自适应进化多目标粒子群优化算法[J]. 控制与决策, 2009, 24(12): 1851-1855 doi: 10.3321/j.issn:1001-0920.2009.12.018Chen Min, Zhang Cong-Yu, Luo Ci-Yong. Adaptive evolutionary multi-objective particle swarm optimization algorithm[J]. Control and Decision Making, 2009, 24 (12): 1851-1855 doi: 10.3321/j.issn:1001-0920.2009.12.018 [35] Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algorithms: NSGA-II. IEEE Trans on Evolutionary Computation, 2002, 6(2):182-197 doi: 10.1109/4235.996017 [36] 熊珞琳, 毛帅, 唐漾, 孟科, 董朝阳, 钱锋. 基于强化学习的综合能源系统管理综述[J]. 自动化学报, 2021, 47(10): 2321-2340Xiong Luo-Lin, Mao Shuai, Tang Yang, Meng Ke, Dong Chao-Yang, Qian Feng. Summary of integrated energy system management based on reinforcement learning[J]. Journal of Automation, 2021, 47 (10): 2321 - 2340 [37] 顾伟, 陆帅, 王珺, 尹香, 张成龙, 王志贺. 多区域综合能源系统热网建模及系统运行优化[J]. 中国电机工程学报, 2017, 5:1305-1315Gu Wei, Lu Shuai, Wang Jun, Yin Xiang, Zhang Cheng-Long, Wang Zhi-He. Heating network modeling and system operation optimization of multi area integrated energy system[J]. Chinese Journal of Electrical Engineering, 2017, 5:1305-1315 [38] 白牧可, 王越, 唐巍, 吴聪, 张博. 基于区间线性规划的区域综合能源系统日前优化调度[J]. 电网技术, 2017, 12:240-247 doi: 10.13335/j.1000-3673.pst.2017.0390Bai Mu-Ke, Wang Yue, Tang Wei, Wu Cong, Zhang Bo. Day ahead optimal scheduling of regional integrated energy system based on interval linear programming[J]. Power Grid Technology, 2017, 12:240-247 doi: 10.13335/j.1000-3673.pst.2017.0390 [39] 吴聪, 唐巍, 白牧可, 张璐, 蔡永翔.基于能源路由器的用户侧能源互联网规划[J]. 电力系统自动化, 2017, 41(04):20-28Wu Cong, Tang Wei, Bai Mu-Ke, Zhang Lu, Cai Yong-Xiang.User side energy Internet planning based on energy router[J]. Power System Automation, 2017, 41 (04): 20-28 期刊类型引用(2)
1. 赵洪山，胡浈，魏伟，温开云. 考虑相关性的配电网分布式光伏承载能力提升方法. 电力系统保护与控制. 2025(01): 37-46 . 
百度学术2. 李浩，李浩，杨扬，朱文超，谢长君. 基于改进鲸鱼算法优化GRU的PEMFC老化预测. 中国电机工程学报. 2024(20): 8166-8178 . 百度学术其他类型引用(2)
-
下载:
计量
- 文章访问数: 629
- HTML全文浏览量: 215
- PDF下载量: 152
- 被引次数: 4
