CN108829946A - 一种基于动态补偿技术的推力计算方法 - Google Patents

一种基于动态补偿技术的推力计算方法 Download PDF

Info

Publication number
CN108829946A
CN108829946A CN201810531926.1A CN201810531926A CN108829946A CN 108829946 A CN108829946 A CN 108829946A CN 201810531926 A CN201810531926 A CN 201810531926A CN 108829946 A CN108829946 A CN 108829946A
Authority
CN
China
Prior art keywords
thrust
dynamic
compensator
cell culture
measuring system
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201810531926.1A
Other languages
English (en)
Inventor
周伟静
常浩
叶继飞
金星
李南雷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Priority to CN201810531926.1A priority Critical patent/CN108829946A/zh
Publication of CN108829946A publication Critical patent/CN108829946A/zh
Pending legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本申请公开了一种基于动态补偿技术的推力计算方法。所述方法通过动态补偿器建立等效系统，利用稳态响应均值计算推力大小。本申请给出了动态补偿器的设计原则，以及基于动态补偿技术的推力计算步骤。该方法不需要修改硬件或进行复杂的硬件调整，适用于所有二阶振动微推力测量系统，具有很强的通用性。

Description

一种基于动态补偿技术的推力计算方法
技术领域
本申请涉及一种基于动态补偿技术的推力计算方法，属于航天器微推进领域的推力测量技术领域。
背景技术
随着新型微推力器的设计、研制和应用，微推力测量技术亦备受关注。微推力测量一般分为稳态推力测量、脉冲平均推力和脉冲冲量测量，稳态推力范围一般在几百微牛至几十毫牛，脉冲冲量范围一般在几十微牛秒至几百微牛秒。稳态推力测量和脉冲冲量测量均采用直接测量方法，即将推力器与推力测量装置固连，利用推力的反冲作用，将推力转化为推力测量装置的力学行为，如振动幅值或转动位移，测量振幅或位移等信息可获得推力或冲量大小。目前，多采用水平旋转的扭摆、吊摆等基于二阶振动系统的微小推力或冲量测量装置。
对于脉冲式微推力器，脉冲冲量测量的实现相对简单，只要力的作用时间小于欠阻尼二阶测量系统振动周期的1/4，可以看做脉冲力瞬间作用于测量系统，根据响应的最大幅值即可计算脉冲冲量。对于稳态微推力器或高频多脉冲微推力器(输出等同于稳态推力)，只要力或多脉冲的作用时间大于欠阻尼二阶测量系统允许误差范围的调整时间，并且维持稳态一定时间，就可以根据响应的稳态幅值计算稳态推力大小。
在机理研究阶段和设计研制过程中，微推力器的推进性能需要测试跟踪，以形成研究闭环，不断改进和提升推进性能。目前，微推力器的研制多处于原理样机或初样机阶段。脉冲式微推力器的单脉冲冲量测量容易实现。但是，由于携带推进剂受限、频率控制技术、功率受限等原因，新型的稳态或多脉冲微推力器，单次试验条件下，开机时长较短，在扭摆响应达到调整时间之前就需关机，无法根据稳态响应还原稳态推力。
通过增大阻尼比和固有频率，可以减小二阶测量系统的调整时间，使响应能尽快进入稳态。但是，一方面，为了能够高精度的标定测量系统的系统参数，阻尼比最好在[0.1,0.4]范围内；另一方面，微推力器的推重比较小，在承载较大质量推力器的基础上，要能够辨识微牛量级的推力，测量系统振动频率均小于1Hz。因此，通过增大阻尼比和固有频率的方法，很难从根本上解决减小调整时间的问题，而且在实验操作时，需要多次调节阻尼比和转动惯量，并多次测试后方可达到测量要求，实施繁琐，周期长。因此，微推力器开机工作时间与测量系统调整时间的匹配问题亟待解决。
发明内容
根据本申请的一个方面，提供了一种基于动态补偿技术的推力计算方法，该方法针对微推力器开机工作时间与测量系统调整时间的匹配难题，提出了补偿器设计原则和补偿实施步骤，可以有效改变测量系统的调整时间，形成新的等效测量系统，使得推力作用在等效测量系统上时，能够利用稳态响应均值计算推力大小。
当微推力器稳态推力持续时间大于二阶机械式直接测量系统周期的0.25倍且小于调整时间时，可通过动态补偿方法建立新的等效测量系统，使得推力作用在等效测量系统上时，能够利用稳态响应均值计算推力大小。本发明给出了动态补偿器的设计原则，以及基于动态补偿技术的推力计算步骤。该方法不需要修改硬件或进行复杂的硬件调整，适用于所有二阶振动微推力测量系统，具有很强的通用性。
所述基于动态补偿技术的推力计算方法，其特征在于，至少包括：
通过动态补偿器建立等效系统，利用稳态响应均值计算推力大小。
可选地，所述动态补偿器的传递函数Hd(s)为：
其中，ωn′为动态补偿器的无阻尼振动频率；
ωn为微推力测量系统的无阻尼振动频率；
s为复变量；
ζ为微推力测量系统的阻尼比；
ζ′为动态补偿器的阻尼比。
可选地，所述动态补偿器满足：
ζ′＝0.7，ωn′＝21.1456/Tf；
其中Tf为稳态力持续时间。
可选地，所述等效系统的建立至少包括：
(a)确定测量装置的基本运动方程；
如采用水平扭摆测量时，水平扭摆是二阶质量-弹簧-阻尼系统，其基本运动方程为方程式(1)；
(b)根据fr＝kθ/b，代入步骤(a)中的基本运动方程，得到测量系统的传递函数；
如将fr＝kθ/b代入水平扭摆的基本运动方程(1)，得到方程式(4)，则获得测量系统的传递函数，即得到方程式(5)；
(c)确认动态补偿器的传递函数；
如在测量系统的传递函数方程式(5)可串连一个阻尼比与固有频率均相同的二阶微分环节，为了减少高频干扰，同样可以加上一个二阶低通滤波器，即将动态补偿滤波器设计成具有方程式(6)所示的传递函数的形式；
(d)获得动态补偿器与测量系统串联后的等效系统；
如方程式(7)所示。通过选取ζ′和ωn′，使得所述等效系统的输出就是符合要求的系统输出。如
ζ′＝0.7；
可选地，所述方法中采用水平扭摆式测量或吊摆测量；
所述动态补偿器为动态补偿滤波器；
所述方法为微小稳态推力计算方法。
可选地，所述方法至少包括：
(1)标定微推力测量系统的系统参数；
(2)在稳态力作用下，采集测量系统的转角输出[ti,θi]，获得平均位置曲线根据fr＝kθ/b，获得
(3)确定动态补偿器参数和传递函数Hd(s)；
(4)将作为Hd(s)的输入，通过设定微推力测量系统的置信区间内端点值，分别获得采样区间内的输出均值Fmax和Fmin，即测量得到的推力为[Fmax,Fmin]；
其中，fr为回复力；
b为推力作用力臂；
k为刚度系数；
θi为时间为ti时的转角。
可选地，步骤(1)中所述微推力测量系统的系统参数包括：
阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k以及置信区间：ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk；
以及测量系统的调整时间ts。
可选地，步骤(3)中所述动态补偿器参数包括：
阻尼比ζ′、无阻尼振动频率ωn′和等效系统的周期
可选地，步骤(3)中所述传递函数Hd(s)为：
其中，ωn′为动态补偿器的无阻尼振动频率；
ωn为微推力测量系统的无阻尼振动频率；
s为复变量；
ζ为微推力测量系统的阻尼比；
ζ′为动态补偿器的阻尼比。
可选地，步骤(4)中所述Fmax的获得方式包括：
将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最大值Fmax；
所述Fmin的获得方式包括：
将作为Hd(s)的输入，当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最小值Fmin。
可选地，步骤(4)中所述采样区间选自[0.8119T',2.4034T']。
可选地，所述方法至少包括：
第一步，标定二阶微推力测量系统的系统参数，获取系统参数阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k及其置信区间ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk，并判断测量系统的调整时间ts；
第二步：推力器的稳态力作用下，采集测量系统的转角输出[ti,θi]，并对数据进行平滑降噪，获得平均位置曲线并根据fr＝kθ/b获得
第三步，确定动态补偿器Hd(s)中的阻尼比ζ′和无阻尼振动频率ωn′，确定传递函数Hd(s)，确定等效系统的周期
第四步，将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最大值Fmax；当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最小值Fmin，则测量得到的推力为[Fmax,Fmin]。
可选地，所述方法至少包括：
首先，利用标定力加载到测量系统上，根据测量系统在标定力下的响应，进行系统动态模型辨识；
然后，根据系统动态模型设计动态补偿滤波器；
最后，将推力器推力加载到测量系统上，将测量系统在推力作用下的响应，经过动态补偿滤波器，还原出推力大小。
可选地，所述方法中的微推力测量系统满足：推力实际工作时间大于微推力测量系统周期的0.25倍，且小于二阶测量系统调整时间。
可选地，所述方法用于基于二阶振动模型的微推力测量系统。
本申请中所述基于动态补偿技术的推力计算方法适用于推力实际工作时间大于二阶测量系统周期的0.25倍，且小于二阶测量系统调整时间的情况。
本申请中所述基于动态补偿技术的推力计算方法适用于所有基于二阶振动模型的微推力测量系统。
作为一种实施方式，本申请中所述方法的具体方案为：
1、基于动态补偿技术的推力计算方法的实施步骤
第一步，标定二阶微推力测量系统的系统参数，获取系统参数阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k及其置信区间ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk，并判断测量系统的调整时间ts；
第二步：推力器的稳态力作用下，采集测量系统的转角输出[ti,θi]，并对数据进行平滑降噪，获得平均位置曲线并根据fr＝kθ/b获得
第三步，确定动态补偿器Hd(s)中的阻尼比ζ′和无阻尼振动频率ωn′，确定传递函数Hd(s)，确定等效系统的周期
第四步，将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最大值Fmax；当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最小值Fmin，则测量得到的推力为[Fmax,Fmin]。
2、动态补偿器的形式及关键参数选取原则
动态补偿器Hd(s)的形式为
动态补偿器Hd(s)中阻尼比ζ′和无阻尼振动频率ωn′的选取原则为：ζ′＝0.7，ωn′＝21.1456/Tf，其中Tf为稳态力持续时间。
3、基于动态补偿技术的推力计算方法的适用情况
基于动态补偿技术的推力计算方法适用于推力实际工作时间大于二阶测量系统周期的0.25倍，且小于二阶测量系统调整时间的情况。
4、基于动态补偿技术的推力计算方法的适用范围
一种基于动态补偿技术的推力计算方法适用于所有基于二阶振动模型的微推力测量系统。
本申请能产生的有益效果包括：
(1)本申请不需要修改硬件或进行复杂的硬件调整，从根本上缩短二阶测量系统的调整时间，大大缩短实验时间。
(2)本申请具有很强的通用性，适用于所有基于水平旋转的扭摆、吊摆等基于二阶振动系统的微小推力测量装置。
附图说明
图1为扭摆系统的典型响应；
图2为本申请一种实施方式中基于动态补偿技术的推力计算步骤。
具体实施方式
下面结合实施例详述本申请，但本申请并不局限于这些实施例。
以水平运动的扭摆测量系统为例，结合附图具体阐述基于动态补偿技术的推力计算方法的提出依据和实施步骤。
1、基于水平扭摆的微小稳态推力测量
扭摆横梁一端安装推力器，另一端安装配重，使得横梁本身、推力器和配重的重心位于转轴上，横梁在水平面内转动。水平扭摆是二阶质量-弹簧-阻尼系统，其基本运动方程为
式中，θ为扭摆横梁相对于横梁零点的转角，J为测量系统转动部件(包括横梁、推力器和配重)相对于转轴的转动惯量，ζ为阻尼比，ωn为无阻尼自振角频率，f(t)为推力器推力，b为推力器推力作用力臂。
因推力器输出的稳态力上升时间非常快，可以当作阶跃力处理，即f(t)＝F，则扭摆系统在阶跃力作用下，系统响应为
式中，当t→∞时，稳态扭转角为从而稳态力仅与系统的稳态扭转角呈线性关系
式中，k为测量系统的扭转刚度系数。
扭摆系统的典型响应如图1所示。ts为扭摆的调整时间，t0为扭摆的1/4振动周期时刻。当稳态力作用时间Tf＝t2，且稳态转角采样区间[ts,t2]大于振动周期的0.2倍，就可以利用采样区间内的稳态转角均值计算稳态推力，且这种平均法带来的相对误差可以控制在一定水平，后面将详细分析。但是当稳态力作用时间Tf＝t1>t0，系统响应还没达到稳态，此时无法利用稳态响应值计算推力。因此，提出利用动态补偿技术改变测量系统的动态特性，减小测量系统的调整时间。
2、动态补偿基本原理及动态补偿滤波器设计
动态补偿的原理是：在原测试系统传递函数上增加一个补偿环节，使得总的传递函数达到理想状态，从而达到改善系统动态特性的目的。补偿环节以辨识建模得到的模型为依据，设计出一种动态补偿滤波器，与原来的测试系统相串联，使级联补偿器后系统的总的动态性能满足使用要求。
根据扭摆测量系统特点，基于动态补偿技术的推力测量方案：首先，利用标定力加载到测量系统上，根据测量系统在标定力下的响应，进行系统动态模型辨识；然后，根据系统动态模型设计动态补偿滤波器；最后，将推力器推力加载到测量系统上，将测量系统在推力作用下的响应，经过动态补偿滤波器，还原出推力大小。
动态补偿滤波器是关键，进行软件补偿即可。常采用的设计方法有直接选择等效系统法和零极点配置法，还有粒子群算法、神经网络等复杂的补偿器设计方法。微推力测量系统属于低阶系统，考虑采用基于直接选择等效系统法的动态补偿滤波器设计方法。
直接选择等效系统法设计动态补偿数字滤波器基本思想是：根据系统动态标定的瞬态响应曲线，由经典控制理论知识，直接得到描述系统特性的动态数学模型，对其模型进行动态分析，当其动态特性不满足测量要求时，给系统串联一个动态补偿数字滤波器，使得串联后构成的等效系统能够满足测量系统的动态特性。
取fr＝kθ/b，则式(1)可改写成
其中，
则测量系统传递函数为
其中，FT＝f(t)，即稳态力；
式(5)可串连一个阻尼比与固有频率均相同的二阶微分环节，为了减少高频干扰，同样可以加上一个二阶低通滤波器，即将动态补偿滤波器设计成具有下述传递函数的形式：
动态补偿滤波器与测量系统串联后的等效系统传递函数为
对于式(5)和式(7)，在系统中加入补偿滤波后，其输出与补偿后等效系统的固有振动频率和阻尼比有关。如果选取合适的ζ′和ωn′，则等效系统的输出就是符合要求的系统输出。下面通过深入分析扭摆稳态响应与稳态推力测量的关系，提出一种ζ′和ωn′的选取方法，提出动态补偿的推力计算步骤。
3、扭摆稳态响应与稳态推力测量的关系
根据式(2)，设在时间区间[T0,T1]内，扭摆系统在阶跃推力作用下的响应均值为
当ζωnT0远大于1且ωd(T1-T0)大于和等于1时，就可以利用稳态转角的均值获取稳态转角，从而利用式(3)计算推力。
利用稳态转角均值计算稳态转角的相对误差为
设且ωd(T1-T0)≥x1≥1，ζωnT0≥x2(且x2远大于1)时，式(9)变为
扭摆系统的高精度标定通常需要阻尼比ζ∈[0.2,0.4]，取阻尼比ζ＝0.3，设g(x1，x2)随x1和x2的变大而迅速变小，当x1≥10和x2≥5时，即k0≥2.5304和k1≥4.122时，g(x1,x2)≤0.0703％。也就是说，在扭摆系统阻尼比ζ＝0.3条件下，如果稳态转角在[2.5304T,4.122T](T为扭摆系统周期)内采样，利用稳态转角均值代替稳态扭转角时，相对误差不超过0.0703％。
根据上述分析，当x1≥10和x2≥5时，扭摆系统阻尼比分别为0.2、0.3、0.4、0.5、0.6、0.7时，稳态转角采样区间分别为[3.8985T,5.49T]、[2.5304T,4.122T]、[1.8233T,3.4149T]、[1.3783T,2.9699T]、[1.061T,2.6526T]、[0.8119T,2.4034T]时，利用稳态转角均值代替稳态扭转角时，相对误差分别≤0.0761％、≤0.0703％、≤0.0682％、≤0.0676％、≤0.0674％、≤0.0674％。因此，具有一定阻尼比的扭摆系统在推力作用下，只要进入稳态的时间区间至少能覆盖某一区间，就可以将稳态转角均值代替稳态扭转角时的相对误差控制在可忽略的水平。
4、动态补偿器设计原则
动态补偿器的设计关键是补偿后等效系统的ζ′和ωn′的选取。根据上述分析，如果稳态力作用时间Tf≤ts2，则可以利用补偿后等效系统的稳态区间采样值计算稳态力。一般，扭摆测量系统阻尼比在[0.2,0.4]范围内，可取补偿后等效系统ζ′＝0.7，则
5、动态补偿下的稳态推力计算的误差来源
由上述分析可知，当推力实际工作时间大于扭摆系统周期的0.25倍，且小于扭摆系统调整时间时，可通过动态补偿方法建立新的等效测量系统，使得推力作用在等效测量系统上时，能够利用稳态响应均值计算推力大小。由上述分析过程可知可知，根据补偿后输出的均值计算稳态力的误差主要来自三个方面：
(1)补偿器输入fr存在的误差
由fr＝kθ/b可知，fr的误差来源于θ、k和b，由于环境等干扰会产生测量噪声，造成实际系统响应数据在平均位置曲线附近上下波动，可采用最小二乘法、正交多项式局部滑动拟合方法等数据处理方法，对实际系统响应测量值进行平滑降噪处理，确定系统响应的平均位置曲线。刚度系数k为具有置信区间的标定值，b存在误差区间的测量力臂值。因此，fr的误差可认为来自于k和b，且k越大、b越小，fr越大。
(2)补偿器Hd(s)中系统参数ζ和ωn带来的误差
系统参数ζ和ωn均为原始测量系统具有置信区间的标定值。其中，ζ主要影响Hd(s)系统稳定的快速性，不影响稳态响应。因此，仅考虑系统参数ωn对补偿器稳态输出的影响，且ωn越大，稳态响应越大。
(3)补偿后输出在[0.8119T′,2.4034T′]区间内取均值带来的误差
在阻尼比为0.7的条件下，在[0.8119T′,2.4034T′]区间内取均值，利用稳态转角均值代替稳态扭转角时的相对误差为≤0.0674％，该误差可以忽略。
6、动态补偿下的稳态推力还原步骤
如图2所示，当推力实际工作时间大于扭摆系统周期的0.25倍，且小于扭摆系统调整时间时，动态补偿下的稳态推力计算步骤如下：
第一步，标定二阶微推力测量系统的系统参数，获取系统参数阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k及其置信区间ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk，并判断测量系统的调整时间ts；
第二步：推力器的稳态力作用下，采集测量系统的转角输出[ti,θi]，并对数据进行平滑降噪，获得平均位置曲线并根据fr＝kθ/b获得
第三步，确定动态补偿器Hd(s)中的阻尼比ζ′和无阻尼振动频率ωn′，ζ′＝0.7，ωn′＝21.1456/Tf，确定传递函数Hd(s)，确定等效系统的周期
第四步，将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最大值Fmax；当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最小值Fmin，则测量得到的推力为[Fmax,Fmin]。
实施例1基于动态补偿技术的推力计算方法
本实施例中所述方法为：通过动态补偿器建立等效系统，利用稳态响应均值计算推力大小。
其中，所述动态补偿器的传递函数Hd(s)为：
其中，ωn′为动态补偿器的无阻尼振动频率；
ωn为微推力测量系统的无阻尼振动频率；
s为复变量；
ζ为微推力测量系统的阻尼比；
ζ′为动态补偿器的阻尼比。
所述动态补偿器满足：
ζ′＝0.7，ωn′＝21.1456/Tf；
其中Tf为稳态力持续时间。
本实施例中采用水平扭摆式测量。
其中，所述等效系统的建立至少包括：
(a)确定测量装置的基本运动方程；
采用水平扭摆测量时，水平扭摆是二阶质量-弹簧-阻尼系统，其基本运动方程为方程式(1)；
(b)根据fr＝kθ/b，代入步骤(a)中的基本运动方程，得到测量系统的传递函数；
将fr＝kθ/b代入水平扭摆的基本运动方程(1)，得到方程式(4)，则获得测量系统的传递函数，即得到方程式(5)；
(c)确认动态补偿器的传递函数；
在测量系统的传递函数方程式(5)可串连一个阻尼比与固有频率均相同的二阶微分环节，为了减少高频干扰，同样可以加上一个二阶低通滤波器，即将动态补偿滤波器设计成具有方程式(6)所示的传递函数的形式；
(d)获得动态补偿器与测量系统串联后的等效系统；
方程式(7)所示。通过选取ζ′和ωn′，使得所述等效系统的输出就是符合要求的系统输出。如
ζ′＝0.7；
其中，利用稳态响应均值计算推力大小包括：
(1)标定微推力测量系统的系统参数；
(2)在稳态力作用下，采集测量系统的转角输出[ti,θi]，获得平均位置曲线根据fr＝kθ/b，获得
(3)确定动态补偿器参数和传递函数Hd(s)；
(4)将作为Hd(s)的输入，通过设定微推力测量系统的置信区间内端点值，分别获得采样区间内的输出均值Fmax和Fmin，即测量得到的推力为[Fmax,Fmin]；
其中，fr为回复力；
b为推力作用力臂；
k为刚度系数；
θi为时间为ti时的转角。
作为一种具体的实施方式，步骤(1)中所述微推力测量系统的系统参数包括：
阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k以及置信区间：ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk；
以及测量系统的调整时间ts。
作为一种具体的实施方式，步骤(3)中所述动态补偿器参数包括：
阻尼比ζ′、无阻尼振动频率ωn′和等效系统的周期
作为一种具体的实施方式，步骤(4)中所述Fmax的获得方式包括：
将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最大值Fmax；
所述Fmin的获得方式包括：
将作为Hd(s)的输入，当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最小值Fmin。
作为一种具体的实施方式，步骤(4)中所述采样区间为[0.8119T',2.4034T']。
实施例2基于动态补偿技术的推力计算方法
本实施例中所述方法与实施例1中区别在于：所述水平扭摆式测量替换为吊摆测量；将水平扭摆测量系统基本运动方程式替换为吊摆测量系统基本运动方程式。
实施例3～7基于动态补偿技术的推力计算方法
实施例3～7中所述方法与实施例1中区别在于：动态补偿器的阻尼比ζ′分别取0.2、0.3、0.4、0.5、0.6，相应的采样区间分别为[3.8985T',5.49T']、[2.5304T',4.122T']、[1.8233T',3.4149T']、[1.3783T',2.9699T']、[1.061T',2.6526T']。
实施例8基于动态补偿技术的推力计算方法
本实施例中所述方法包括：
第一步，标定二阶微推力测量系统的系统参数，获取系统参数阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k及其置信区间ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk，并判断测量系统的调整时间ts；
第二步：推力器的稳态力作用下，采集测量系统的转角输出[ti,θi]，并对数据进行平滑降噪，获得平均位置曲线并根据fr＝kθ/b获得
第三步，确定动态补偿器Hd(s)中的阻尼比ζ′和无阻尼振动频率ωn′，确定传递函数Hd(s)，确定等效系统的周期
第四步，将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最大值Fmax；当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最小值Fmin，则测量得到的推力为[Fmax,Fmin]。
其中，所述动态补偿器Hd(s)的形式为
所述动态补偿器Hd(s)中阻尼比ζ′和无阻尼振动频率ωn′的选取原则为：ζ′＝0.7，ωn′＝21.1456/Tf，其中Tf为稳态力持续时间。
实施例9基于动态补偿技术的推力计算方法
本实施例中所述方法包括：
首先，利用标定力加载到测量系统上，根据测量系统在标定力下的响应，进行系统动态模型辨识；
然后，根据系统动态模型设计动态补偿滤波器；
最后，将推力器推力加载到测量系统上，将测量系统在推力作用下的响应，经过动态补偿滤波器，还原出推力大小。
以上所述，仅是本申请的几个实施例，并非对本申请做任何形式的限制，虽然本申请以较佳实施例揭示如上，然而并非用以限制本申请，任何熟悉本专业的技术人员，在不脱离本申请技术方案的范围内，利用上述揭示的技术内容做出些许的变动或修饰均等同于等效实施案例，均属于技术方案范围内。

Claims (10)

1.一种基于动态补偿技术的推力计算方法，其特征在于，至少包括：
通过动态补偿器建立等效系统，利用稳态响应均值计算推力大小。
2.根据权利要求1所述的方法，其特征在于，所述动态补偿器的传递函数Hd(s)为：
其中，ω′n为动态补偿器的无阻尼振动频率；
ωn为微推力测量系统的无阻尼振动频率；
s为复变量；
ζ为微推力测量系统的阻尼比；
ζ′为动态补偿器的阻尼比；
优选地，所述动态补偿器满足：
ζ′＝0.7，ω′n＝21.1456/Tf；
其中Tf为稳态力持续时间。
3.根据权利要求1所述的方法，其特征在于，所述方法中采用水平扭摆式测量或吊摆测量；
所述动态补偿器为动态补偿滤波器；
所述方法为微小稳态推力计算方法。
4.根据权利要求1所述的方法，其特征在于，所述方法至少包括：
(1)标定微推力测量系统的系统参数；
(2)在稳态力作用下，采集测量系统的转角输出[ti,θi]，获得平均位置曲线根据fr＝kθ/b，获得
(3)确定动态补偿器参数和传递函数Hd(s)；
(4)将作为Hd(s)的输入，通过设定微推力测量系统的置信区间内端点值，分别获得采样区间内的输出均值Fmax和Fmin，即测量得到的推力为[Fmax,Fmin]；
其中，fr为回复力；
b为推力作用力臂；
k为刚度系数；
θi为时间为ti时的转角。
5.根据权利要求4所述的方法，其特征在于，步骤(1)中所述微推力测量系统的系统参数包括：
阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k以及置信区间：ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk；
以及测量系统的调整时间ts；
优选地，步骤(3)中所述动态补偿器参数包括：
阻尼比ζ′、无阻尼振动频率ω′n和等效系统的周期
优选地，步骤(3)中所述传递函数Hd(s)为：
其中，ω′n为动态补偿器的无阻尼振动频率；
ωn为微推力测量系统的无阻尼振动频率；
s为复变量；
ζ为微推力测量系统的阻尼比；
ζ′为动态补偿器的阻尼比。
6.根据权利要求4所述的方法，其特征在于，步骤(4)中所述Fmax的获得方式包括：
将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最大值Fmax；
所述Fmin的获得方式包括：
将作为Hd(s)的输入，当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在采样区间内取输出均值即测量得到的推力最小值Fmin；
优选地，步骤(4)中所述采样区间选自[0.8119T',2.4034T']。
7.根据权利要求4所述的方法，其特征在于，所述方法至少包括：
第一步，标定二阶微推力测量系统的系统参数，获取系统参数阻尼比ζ、无阻尼振动频率ωn、周期T和刚度系数k及其置信区间ζ±Δζ、ωn±Δωn、T±ΔT、k±Δk，并判断测量系统的调整时间ts；
第二步：推力器的稳态力作用下，采集测量系统的转角输出[ti,θi]，并对数据进行平滑降噪，获得平均位置曲线并根据fr＝kθ/b获得
第三步，确定动态补偿器Hd(s)中的阻尼比ζ′和无阻尼振动频率ω′n，确定传递函数Hd(s)，确定等效系统的周期
第四步，将作为Hd(s)的输入，当k＝k+Δk、b＝b-Δb、ωn＝ωn+Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最大值Fmax；当k＝k-Δk、b＝b+Δb、ωn＝ωn-Δωn时，获得补偿器的输出，从推力作用时刻开始，在[0.8119T′,2.4034T′]区间内取输出均值即测量得到的推力最小值Fmin，则测量得到的推力为[Fmax,Fmin]。
8.根据权利要求1所述的方法，其特征在于，所述方法至少包括：
首先，利用标定力加载到测量系统上，根据测量系统在标定力下的响应，进行系统动态模型辨识；
然后，根据系统动态模型设计动态补偿滤波器；
最后，将推力器推力加载到测量系统上，将测量系统在推力作用下的响应，经过动态补偿滤波器，还原出推力大小。
9.根据权利要求1所述的方法，其特征在于，所述方法中的被测推力满足：推力实际工作时间大于微推力测量系统周期的0.25倍，且小于二阶测量系统调整时间。
10.根据权利要求1所述的方法，其特征在于，所述方法用于基于二阶振动模型的微推力测量系统。
CN201810531926.1A 2018-05-29 2018-05-29 一种基于动态补偿技术的推力计算方法 Pending CN108829946A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810531926.1A CN108829946A (zh) 2018-05-29 2018-05-29 一种基于动态补偿技术的推力计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810531926.1A CN108829946A (zh) 2018-05-29 2018-05-29 一种基于动态补偿技术的推力计算方法

Publications (1)

Publication Number Publication Date
CN108829946A true CN108829946A (zh) 2018-11-16

Family

ID=64146649

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810531926.1A Pending CN108829946A (zh) 2018-05-29 2018-05-29 一种基于动态补偿技术的推力计算方法

Country Status (1)

Country Link
CN (1) CN108829946A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109579969A (zh) * 2018-11-29 2019-04-05 上海交通大学 叶轮在加减速瞬态工况下最大振动幅值的获取方法及系统
CN109682526A (zh) * 2019-03-04 2019-04-26 合肥工业大学 一种多维力传感器的动态解耦方法
CN110413015A (zh) * 2019-06-27 2019-11-05 北京控制工程研究所 基于闭环控制的微牛量级微推力动态测试台及测试方法
CN112231890A (zh) * 2020-09-03 2021-01-15 兰州空间技术物理研究所 基于扭摆测量系统的高平稳电推力器的推力测评方法
CN114218750A (zh) * 2021-11-16 2022-03-22 中国人民解放军战略支援部队航天工程大学 基于数字滤波器的星载微推力器推力响应时间测量方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104990666A (zh) * 2015-05-27 2015-10-21 中国人民解放军装备学院 一种基于比例回归法的二阶振动测量系统的系统参数标定方法
CN105784237A (zh) * 2016-05-13 2016-07-20 中国科学院力学研究所 一种微推力测试系统及方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104990666A (zh) * 2015-05-27 2015-10-21 中国人民解放军装备学院 一种基于比例回归法的二阶振动测量系统的系统参数标定方法
CN105784237A (zh) * 2016-05-13 2016-07-20 中国科学院力学研究所 一种微推力测试系统及方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
周伟静 等: ""一种基于动态补偿技术的微小稳态推力还原方法"", 《航空学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109579969A (zh) * 2018-11-29 2019-04-05 上海交通大学 叶轮在加减速瞬态工况下最大振动幅值的获取方法及系统
CN109579969B (zh) * 2018-11-29 2020-05-22 上海交通大学 叶轮在加减速瞬态工况下最大振动幅值的获取方法及系统
CN109682526A (zh) * 2019-03-04 2019-04-26 合肥工业大学 一种多维力传感器的动态解耦方法
CN110413015A (zh) * 2019-06-27 2019-11-05 北京控制工程研究所 基于闭环控制的微牛量级微推力动态测试台及测试方法
CN112231890A (zh) * 2020-09-03 2021-01-15 兰州空间技术物理研究所 基于扭摆测量系统的高平稳电推力器的推力测评方法
CN114218750A (zh) * 2021-11-16 2022-03-22 中国人民解放军战略支援部队航天工程大学 基于数字滤波器的星载微推力器推力响应时间测量方法

Similar Documents

Publication Publication Date Title
CN108829946A (zh) 一种基于动态补偿技术的推力计算方法
CN108897226B (zh) 基于干扰观测器的mems陀螺仪预设性能非奇异滑模控制方法
CN102914972B (zh) 基于模型整体逼近的微陀螺仪rbf网络自适应控制方法
CN102426420B (zh) 一种强鲁棒性的运动载体光电稳定平台控制系统
KR101223669B1 (ko) 엔진 벤치 시스템 제어 시스템
Park et al. Adaptive control for the conventional mode of operation of MEMS gyroscopes
CN103116275B (zh) 基于滑模补偿的微陀螺仪鲁棒神经网络控制系统及方法
CN107607101B (zh) 基于干扰观测器的mems陀螺滑模控制方法
CN105043348A (zh) 基于卡尔曼滤波的加速度计陀螺仪水平角度测量方法
CN105159084B (zh) 一种带干扰观测器的机械手神经网络控制系统及方法
CN103324087B (zh) 基于神经网络的微陀螺仪的自适应反演控制系统及方法
RU2285902C1 (ru) Способ определения и компенсации ухода гиростабилизированной платформы и устройство для его осуществления
CN108205259A (zh) 基于线性扩张状态观测器的复合控制系统及其设计方法
CN115630558B (zh) 一种复合材料构件装配变形预测方法
CN103345155B (zh) 微陀螺仪的自适应反演控制系统及方法
CN110389528A (zh) 基于扰动观测的数据驱动mems陀螺仪驱动控制方法
CN106842953B (zh) 一种无人直升机自适应低阶控制器
CN107844618A (zh) 用于测量推力和冲量的扭摆系统的设计方法
NO339405B1 (no) Fremgangsmåte og anordning for å styre/regulere en fysisk variabel i et dynamisk system.
CN108108523A (zh) 一种飞机荷兰滚等效拟配初值选取方法
CN109062048A (zh) 基于复合学习的mems陀螺仪预设性能非奇异滑模控制方法
EP3250886B1 (en) Gyroscope loop filter
CN109426143A (zh) 负载转矩估算方法、系统、机电控制系统、方法及电机
CN112799303B (zh) 一种机械臂的h∞控制方法
CN107861384B (zh) 基于复合学习的mems陀螺仪快速启动方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20181116