复杂系统的设计已经证实了仿真是从概念设计到样品实现过程中不可或缺的步骤。事物仿真需要考虑初步估计、择优甚至是现有方案在实施前的重新设计,因此会降低风险。建立元器件本身及其之间交互作用的高保真模型,能够使仿真的结果更具有确信度。

联合仿真是现阶段最有效的耦合系统仿真方案。从数学方面分析,联合仿真由一组常微分方程和微分代数方程构成,而仿真时间包括子系统之间的通信节点时间和每一个子系统内部的积分时间。分析一个联合仿真系统是否符合要求需要对其稳定性、准确性和快速性进行分析。

联合仿真系统中,每一个仿真器都可以看成是一个动态系统,它们之间具有信号的传递,传递的信号与上一次信号传递的时间间隔为一个宏步长,该步长的大小反映了系统仿真计算时间的长短。通常情况下,仿真过程结果都希望具有最小的通信时间,然而考虑到系统稳定性,该计算时间不能无限变小,所以需要找出兼顾计算效率和计算稳定性的方法来满足系统在仿真需求。为了提高仿真速度,联合仿真的并行化速度也需要提高,然而并行的子模型和相关求解器之间的松耦合会引起积分误差,为了保证误差在预先要求的范围内,同样需要兼顾考虑仿真效率和仿真精确度的关系。

Dymola是一款基于Modelica语言的多领域仿真平台。Modelica语言可以实现各模块库之间的连接,进行多领域仿真。文献[1]通过Modelica语言搭建了包括交流发电机、变压整流器、电源系统接触器和过欠压、过欠频检测模块等部件的数字仿真模型和供电系统模型,并完成系统逻辑设计。通过特殊的接口,Dymola可以与多款软件实现联合仿真,现在国内外的主要研究方向为对电机的控制和汽车系统建模。文献[2]利用Simulink、CarSim、Dymola软件实现对负荷传感液压转向系统的联合仿真;Dymola符合FMI接口协议,所以可以输入输出FMU模型;文献[3]在Dymola中实现二极管模型的FMU模型输出;文献[4]在Dymola中实现热流体模型的FMU模型输出;文献[5]是对自动嵌入式物理模型的FMI研究;文献[6]将Modelica模型通过FMI输出到Excel表格中,通过FMI Toolbox for Matlab工具生成FMU模型,最后导入Dymola平台中进行联合仿真,实现控制器算法模型和Dymola整车联合仿真测试。

本文对联合仿真系统模型进行了数学分析,建立了系统微分方程,在稳定和精度要求下提出了线性隐性计算方法,并提出了基于FMI协议的联合仿真计算流程,最后通过一个工程例子说明了该方法的实际应用效果。

用非线性微分方程描述一个连续离散混合动态系统Σ形式如下

式中,$\text{X}\in {{R}^{{{n}_{x}}}}$是连续状态向量;$\text{D}\in {{R}^{{{n}_{D}}}}$是离散状态向量;$\text{U}\in {{R}^{U}}$是输入向量;$\text{Y}\in {{R}^{{{n}_{Y}}}}$是输出向量;$t\in {{R}^{+}}$表示时间。

严格增加的时间变量序列${{({{t}_{n}})}_{n\ge 0}}$代表被称为状态事件的间断点,它们是$h\left( t,X,D,U \right)=0$方程的根。函数h通常被称为事件指示器或过零函数,用于进行事件检测和定位。

每一个时刻都会计算得到一个连续状态向量$X\left( {{t}_{n}} \right)=I\left( {{t}_{n}},X,D,U \right)$作为事件句柄的结果和一个离散状态变量$D\left( {{t}_{n}} \right)=J\left( {{t}_{n-1}},X,D,U \right)$作为离散状态更新的结果。如果$X\left( {{t}_{n}} \right)$没有受到间断点的影响,该向量的右极限将会和它在${{t}_{n}}$处的取值相等。

假设Σ对于变量$X\left( {{t}_{0}} \right)$和$D\left( {{t}_{0}} \right)$均存在相应的初始状态,从而使得X、D、U和Y在$\left[ {{t}_{n}},{{t}_{n+1}} \right]$内是分段光滑的连续函数。将模型分为几个子模型使得系统能够并行运行,为了简化分解过程,图1将整个系统分解成了两个子系统。该方法可以推广至将系统Σ任意分解为N个模块^[7]。

因此系统Σ可以重新写成

式中,${{U}_{1}}$是子系统${{\sum }_{1}}$的输入,${{U}_{2}}$是子系统${{\sum }_{2}}$的输入。也就是说,${{U}_{1}}\bigcup {{U}_{2}}=U$,${{U}_{1}}\bigcap {{U}_{2}}$根据得到的解耦结果判断是不是空集。同样,${{Y}_{1}}$是子系统${{\sum }_{1}}$得到的输出,${{Y}_{2}}$是子系统${{\sum }_{2}}$得到的输出,即${{Y}_{1}}\bigcup {{Y}_{2}}=Y$,${{Y}_{1}}\bigcap {{Y}_{2}}=\varnothing $。

为了更好数值运算整个复杂系统,每一个仿真器都需要在通信节点进行数据交换,来为其他子模块提供所需要的数据。为了提高积分运算速度,并行模块要尽量独立,所以子模型之间的同步时间P要远大于他们内部的积分时间${{h}_{n}}\left( P\gg {{h}_{n}} \right)$。因此在通信节点之间,每一个仿真器都有自己的积分速率(假设使用变步长求解器),并且在此期间,从其他模型处得到的数据将保持常数。

大的通信间隔很可能提高积分运算速度,但是也可能会产生积分误差累计从而降低最终结果的精确度。所以通过非严格同步来减小建模误差是找到提高积分速度和精确度之间平衡的有效方式的第 一步。

对于图2中的模型,每一个子系统的输出变量在通信节点${{t}_{s}}\left( {{t}_{s+1}}-{{t}_{s}}={{H}_{s}} \right)$处被采样,该值并作为输入变量在整个宏步长${{H}_{s}}$期间保持常数,所以整体看来,每一个子系统中的求解器是一个离散动态系统。对于宏步长${{H}_{k}}$的频率选取,根据香农采样定理,需要使其至少大于Nyquist频率的二倍,即采样角频率${{\omega }_{samp}}$满足

式中,${{\omega }_{\max }}$为系统最大角频率,满足

所以系统的宏步长满足

然而这个宏步长的边界太大,由经验求得,宏步长的频率要远大于Nyquist频率,有时会达到100倍。如果宏步长的频率不够大,则会导致系统的最终发散。为了解释只有很小的宏步长才能使耦合系统稳定,需要对整个系统的稳定性进行分析,按照文献[8]中所述的方法可以对这类环路采样系统进行稳定性分析。此时,子系统通过零阶保持器来耦合变量。

图3为两个子系统的联合仿真方框图,其主要在每一步离散化了耦合变量,这模拟了子系统可以分别使用数值求解器进行内部计算。对于两个系统的耦合连接,通过线性化状态空间所求出的离散时间传递函数来分析稳定性^[9]。

若要保证系统稳定,需要足够小的宏步长,这虽然保证了系统的精确度,但却需要消耗过多的计算时间,降低了计算效率^[10]。线性隐性稳定方法利用子系统的Jacobian矩阵来线性化系统,这种方法允许使用相对较大的宏步长来实现系统的稳定^[11]。

整个耦合系统的数学模型如下所示

式中,$x\left( t \right)=\left( {{x}^{1}}\left( t \right),{{x}^{2}}\left( t \right),\cdots,{{x}^{N}}\left( t \right) \right)$表示N个子系统的状态变量集合,$u\left( t \right)=\left( {{u}^{1}}\left( t \right),{{u}^{2}}\left( t \right),\cdots,{{u}^{N}}\left( t \right) \right)$和$y\left( t \right)=\left( {{y}^{1}}\left( t \right),{{y}^{2}}\left( t \right),\cdots,{{y}^{N}}\left( t \right) \right)$分别表示子系统输入输出集合。

对于整个系统,内部子系统之间的连接是通过耦合方程来实现的。

在耦合方程中,每一个子系统的输入向量是其他子系统输出向量的函数,K为具有以下特点的关联矩阵：①K是一个方阵,认为每一个子系统的输出都对应一个子系统的输入;② 只包含1或者0;③ 矩阵中每一行和每一列中只有一个1。

对系统进行线性化处理,规定Jacobian矩阵为

在使用隐性稳定方法时,需要以下假设：① 在没有考虑耦合变量的代数环情况下,假设DK的乘积为幂零;② 在联合仿真宏步长$\left( t\in \left( {{t}_{k}},{{t}_{k}}+{{H}_{k}} \right) \right)$中,线性时不变系统由线性化(在点$x=x\left( {{t}_{k}} \right)$和$u=u\left( {{t}_{k}} \right)$系统常微分方程得到;③ 通过Backward Euler方法来离散线性化后的系统。

根据假设②,线性化后的系统数学模型如下 所示

式中,$\xi $、$\eta $和$\omega $分别在线性化系统中对应于x、y和u的值。因此,在一个宏步长内,重新描述线性化处理后的系统耦合方程为

式(9)为整个系统的微分代数方程,因为不考虑代数环,微分代数方程可以通过下面的计算改写成常微分方程。根据假设③,常微分方程首先进行离散化处理

Backward Euler:

$\dot{\xi }\left( t \right)\approx {{\left( I-\left( t-{{t}_{k}} \right)A \right)}^{-1}}\left[ \dot{x}\left( {{t}_{k}} \right)+ B\left( \omega \left( t \right)-u\left( {{t}_{k}} \right) \right) \right]$

同样可以得到输出的关系式

Backward Euler:

$\begin{matrix} \dot{\eta }\left( t \right)\approx C\dot{\xi }\left( t \right)+D\dot{\omega }\left( t \right) \\ \dot{\eta }\left( t \right)\approx C{{\left( I-\left( t-{{t}_{k}} \right)A \right)}^{-1}}\left[ \dot{x}\left( {{t}_{k}} \right) \right.\left. B\left( \omega \left( t \right)-u\left( {{t}_{k}} \right) \right)+D\dot{\omega }\left( t \right) \right] \\ \end{matrix}$

为了考虑缺少代数环,式(10)可以通过假设②估计得到,即输入ω是通过在$t={{t}_{k}}$处线性化系统的输出y得到的。

Backward Euler:

$\dot{\omega }\left( t \right)=K\left[ -C\dot{x}\left( t \right)+DK\dot{\eta }\left( t \right) \right]$

将式(13)代入到式(12)中,并根据假设①设定DK的乘积为幂零,得

故式(10)可以写成

根据式(16),每一步结束后输出值会在子系统中进行求解和传递

每个功能模型单元(Function mockup unit,FMU)只负责本地系统运算,不包含任何耦合和仿真环境的相关信息^[12]。协同仿真过程中,计算机只通过主求解器来联合系统中不同从求解器的信息,并向从求解器提供方法来求解其他代表耦合系统的线性部分的常微分方程(Ordinary differential equation,ODE)。这些信息一部分来自于每一个求解器的模型描述结构中,一部分来自于求解FMU所包含的系统方向导数的函数^[13]。

FMU求解器需要具有计算所有状态变量和所连的输入输出的方向导数的能力,在FMI协议中,通过标记“provide Direction Derivative”来说明,并通过方程“fmi Get Directional Derivative”来实现。稳定方法“基于模型的外推法”不与FMI提供的“基于历史的外推法”相兼容,所以属性“can Interpolate Inputs”需要扩展到其所指定的外推类型。

将DOE转变为ODE的计算流程为首先建立包含s$\left( s\in \left\{ 1,2,\cdots,N \right\} \right)$个从求解器的矩阵方程：定义方阵${{K}_{\bar{s},s}}$为纵坐标作为当前求解器的输出变量,横坐标为来自其他求解器的输入变量;定义方阵${{K}_{\bar{s},s}}$为横坐标作为当前求解器的输入变量,纵坐标为来自其他求解器的输出变量。

计算流程包括两步,第一步在通信节点处进行,第二步发生在求解在通信节点之间的连续时间DOE系统处。

第一步中状态向量的倒数求得为

它的输出为

Jacobian矩阵为

在临界值$t_{k}^{+}$处的值为

上式中当k等于0时,表示为全局的初始状态。

在${{t}_{k}}<t\le {{t}_{k+1}}$时,求解耦合常微分方程系统的一个子集为

其中,${{f}_{{\bar{s}}}}$由主求解器求得

根据式(19)可知,出现在式(21)中的输入变量包含了其他的状态变量,并且输入变量现在相当是状态变量的初始状态。通过求解器内的数值积分求解出的状态向量是$\left( {{x}_{s}},{{u}_{s}} \right)$,这里${{u}_{s}}$实际是${{K}_{s,\bar{s}}}{{\eta }_{{\bar{s}}}}$,通过主求解器来连接线性系统的输出${{\eta }_{{\bar{s}}}}$和求解器的输入,因此连接的结构信息只用主求解器获得。

如果主求解器不支持稳定性,式(22)的右边会降为0,这样从求解器的实际输入将不变：如果没有设置标记“can Interpolate Inputs”,那么零阶保持器外推因主求解器不支持稳定性而在FMI协议中定义。这种情况下,通过主求解器在每一个通信节点$\left( {{t}_{k}},k=0,\cdots,M-1 \right)$调用方程“fmi Get Real”和“fmi Set Real”,每一个从求解器的输入是通过式(20)直接给出的,${{y}_{{\bar{s}}}}$是从求解器的输出,这些从求解器会求出第s个求解器的输入。

主求解器方面需要设置回调方程(Callback function),这个方程用于求解相连的从求解器的动态输入。在联合仿真之前,从求解器的模型结构信息将体现在准备矩阵和式(23)的计算中。随后,在联合仿真过程中,这些元素将会调用方程“fmi Get Continuous States”“fmi Get Derivatives”和“fmi Get Directional Derivatives”在每一次的信息交换节点上更新。当对式(21)、(23)数值积分时,系统通过从数值求解器调用这些方程。

根据式(23)可知,回调方程需要返回稳定输入的微分向量,这需要通过主仿真环境在从设备实现方程属性“fmi Callback Functions”来实现。图4简略地描述了主求解器中的算法流程^[14]。

利用FMI2.0 for Co-simulation标准的线性隐性稳定方法可以很好对系统加速,其改善了联合仿真在平衡稳定或精确度和计算效率上的问题。随着FMI2.0的问世,一个巨大的提升是出现了方向导数矩阵的接口,这使得联合仿真系统的强耦合性更有效^[15,16]。

基于二自由度液压模型如图5所示,其中包括了两个液压子系统,可将这个系统看为两个子系统实现联合仿真。

所有测试都采用同一套参数,测试结果如表1所示。

联合仿真的精度是通过方均根误差来表示的。在所有的联合仿真测试中,CPU时间通过设定一定的间隔时间为50 μs。从表1中可以看出,混合使用C/Python标准的#5方案相比#3方案在CPU耗时上更差,这是由于在两个仿真环境中带有很多的转换接口;相反,对比#2和#6,同样使用C标准的#3和#7(或#8),系统明显实现了加速。说明利用FMI2.0标准的线性隐性稳定方法可以很好地改善联合仿真在平衡稳定或精确度和计算效率上的问题。

本文的研究目的是提高联合仿真系统数值积分运算速度,同时兼顾系统稳定运行并保证积分误差在预期范围内。论文研究了一种基于FMI通用接口标准的联合仿真方法,得到如下结论。

(1) 其基本方法是将整个系统分为若干个子系统,各部分之间在通信节点处交换数据,数值分析时需要将系统合理解耦分析,在保证运算速度的同时使系统具有良好的稳定性和准确性。

(2) 线性隐性稳定方法在保证系统稳定的前提下对仿真运算时间具有一定的加速,其基于联合仿真的FMI标准下的工程实例也证明了该方法使系统具有良好的稳定性和准确性。