复杂产品多学科模型的协同仿真求解算法研究

第21卷第20期系

统仿真学报? V ol. 21 No. 20

2009年10月 Journal of System Simulation Oct., 2009

梁思率，张和明

(清华大学国家CIMS 工程技术研究中心，北京 100084)

摘要：在数值计算领域求解大型常微分方程组时通常采用分解算法，类似地，在多学科协同仿真中系统模型往往被拆分成多个子模型并采用多个求解器进行求解。基于以上相关性，研究了多学科协同仿真算法的基本原理，在微分方程组合算法的基础上提出了基于联合仿真步的组合算法，给出了算法的形式化描述和原理说明，并通过一个具体实例验证了算法的有效性。关键词：虚拟样机；多学科协同仿真；微分方程；组合算法

中图分类号：TP391.9 文献标识码：A 文章编号：1004-731X (2009) 20-6368-05

Research on Algorithm of Collaborative Simulation

for Multidisciplinary Model of Complex Product

LIANG Si-lv, ZHANG He-ming

(National CIMS Engineering Research Center, Tsinghua University, Beijing 100084, China)

Abstract: It is common to use decomposition methods to solve large ordinary differential equations in numerical calculation. Similarly, system model is often divided into several sub-models and solved by several solvers in multidisciplinary collaborative simulation. Based on this correlation, basic idea of algorithm for multidisciplinary collaborative simulation was studied. Combinative algorithm based on united simulation step was proposed based on combinative algorithm for ordinary differential equations. Formalized description and theoretical illustration of the algorithm were given. Finally an example was used to check validity of the algorithm.

Key words: virtual prototyping; multidisciplinary collaborative simulation; ordinary differential equations; combinative algorithm

引言复杂产品已发展成为以机械系统为主体，融合电子、控制、液压、气动和软件等多领域集成的物理系统，机电液控等多领域耦合是当前复杂产品虚拟样机系统的重要特点。协同仿真逐渐成为多领域耦合复杂机械产品开发的重要环节，是提高其动态特性与综合性能的关键技术。复杂产品的仿真模型往往是由多个子模型组成的，每个子模型属于一个单独的学科(动力学、控制、液压等)。多学科协同仿真即指通过多个属于不同学科的子模型之间的并行运行、实时交互，共同实现对涉及学科间交互的复杂产品设计问题的仿真分析

[1]

。通常有两种方法：一是基于各种单学科建模仿真软件，

通过软件之间的接口或仿真总线，来完成虚拟样机的多学科建模和协同仿真；二是针对物理系统具体特点，通过数学抽象的方式实现多学科系统的统一描述和模型的求解。后者的典型例子是Modelica 语言，由于不涉及系统分解，因此是对完整模型采用单一算法进行仿真求解，计算精度高，但是多学科模型的统一建模比较不易，且没有利用已有的各学科专用建模软件的功能。前者目前有两种常见的实现方式：基

收稿日期：2008-06-16 修回日期：2008-08-14

基金项目：国家自然科学基金 (60674079)；国防基础科研项目 (B0420060524) 作者简介：梁思率(1984-), 男, 浙江杭州人, 硕士, 研究方向为虚拟样机, 多学科协同仿真；张和明(1966-), 男, 浙江仙居人, 副教授, 研究方向为虚拟样机, 并行工程。

于联合仿真接口的协同仿真[2,3]和基于RTI 的分布式协同仿真[4,5]。无论采用哪一种，其仿真原理都是将每个子模型交由各自的仿真软件进行解算，子模型之间则采用一定的交互机制完成输入输出数据的更新。由于模型通常为连续系统，即可表示成微分方程组的形式，因此仿真的本质是一数值积分过程，相应的仿真软件称为积分器。一个完整的模型分解成若干个子模型后，原先采用单个积分器求解的过程就转化为采用多个积分器协同求解的过程。在数值计算领域求解大型常微分方程组时往往采用分解算法，即将其分拆成多个较小的微分方程组分别进行求解，其算法的收敛性是得到证明的[6]。多学科协同仿真采用的是类似的思想，本文研究了多学科协同仿真算法的基本原理，及其模型求解时的组合算法，并通过具体实例验证了该算法的有效性。

1 协同仿真算法原理

1.1模型描述

连续系统仿真中的模型一般可以描述成微分方程组的形式。

假设一个完整系统的模型可以描述为大型常微分方程

组的初值问题，即：

00()

((),),(),M N dz t f z t t z t z z R dt +==∈ (1) 其中，t 是时间，()z t 是状态向量。对模型的仿真过程就是对该微分方程组的求解过程。

2009年10月梁思率，等：复杂产品多学科模型的协同仿真求解算法研究 Oct., 2009

在复杂产品的多学科协同仿真建模过程中，当一个整体模型(系统)被分解成多个子模型(子系统)后，相当于一个完整的微分方程组被拆分成多个相互耦合的子微分方程组。为简单起见，考虑将其分解成两个子模型的情形：

111211010112112112()((),(),),(),()((),(),),M

S dz t f z t Y t t z t z z R dt

Y t y z t Y t t Y R

?==∈?

??=∈?

222122020221

221221()

((),(),),(),()((),(),),N T

dz t f z t Y t t z t z z R dt

Y t y z t Y t t Y R ?==∈?

??=∈? (3) 其中，t 是时间，12(),()z t z t 分别是子模型1、2的状态向量，

1221(),()Y t Y t 分别是子模型1、2的输出向量，同时也分别是子模型2、1的输入向量。

1.2 微分方程的组合算法[6]

对完整微分方程组的数值积分求解过程就转化为对多个子微分方程组运用不同的数值积分器进行求解的过程，类似于数值计算中采用的基于系统分割的组合算法[6]。

考虑一个常微分方程组描述的系统：

00(,),(),M N dz

f z t z t z z R dt

+==∈将其分解成两个子系统S 1、S 2：

111210101(,,),(),M z

f z z t z t z z R ==∈ 222120202(,,),(),N z

f z z t z t z z R ==∈ 假设子系统S 1是快变系统，其积分方法为F ，积分步长为h ；子系统S 2是慢变系统，其积分方法为S ，积分步长为,,0H H rh r =>。在计算一个子系统的右端函数时耦合变量用另一个子系统输出结果的插值代替，分别表示为

12,z

z ，其中，子系统S 1采用插值公式I ，子系统S 2采用插值公式J ，插值公式可以是Lagrange 多项式或Hermite 多项式。假设,m F t 是方法F 的结点，,n s t 是方法S 的结点，则在,n s t 和1,n s t +之间有方法F 的1r ?个结点，如图1所示。

图1 多速率系统积分过程

假设,n s t 时刻方法S 和F 的数值解分别为2,2,1,,,,n s n s rm F z f z 1,rm F f ，其中，,1,(1),,n s n s r m F rm F t t H t rh t ??=+=+=。一个通用

的组合算法如图2所示。

该算法实际上是一种并行算法，因所有第1i +个循环内的数值解均得自第i 个循环及以前的解的插值。因此，步骤

2（第2个框）到步骤6和步骤7到步骤8可以同时执行。图2 通用组合算法

1.3 组合算法的收敛性

费景高[7]证明了该组合算法能收敛到系统分解前微分方程组的真实解，并研究了其收敛阶，有如下定理：

定理1：如果积分方法F ，S 和插值公式I ，J 的阶分别是p 1，p 2，p 3，p 4的话，则算法的收敛阶为p ，p 是它们之中最小者。

该定理表明，插值公式阶数的选取要根据积分方法的阶数来确定，若插值公式的阶数取得太高，则容易引起积分不稳定性，若取得太低，则会使算法的收敛阶降低。

该定理对其适用范围作了如下限定：系统方程采用微分方程描述，不能出现分段函数等不连续的情况，而实际建模过程中也应尽量避免使用非连续函数对系统进行描述[8]；系

统分解的方式可以是刚性也可以是非刚性的，分解的个数可以是2个或者更多；积分的方法可以是单步法也可以是多步法；该定理虽然是基于等步长积分方法而作出的，但同样适用于变步长方法的场合。这一结论可以有效支持类似基于这种分解方式的组合算法，包括复杂产品多学科协同仿真的求解过程。

然而，对于复杂产品虚拟样机的多学科协同仿真而言，首先，其多学科系统的分解方式实际上并非是基于快慢子系

S 2

S 1

t n,S (t rm ,F ) t n +1,S (t r (m+1),F ) S,J

F ,I

2009年10月系统仿真学报 Oct., 2009

统的划分，它们通常事先分解成动力学、控制、气动等多个子系统，相互之间并不存在快慢关系，采用的积分方法也是多种多样，比如，可以是变阶变步长的(如BDF 算法)；其次，定理1的证明过程是基于常微分方程(ODE)表示的模型，而实际的产品模型可以是基于微分代数方程(DAE)的描述，其积分的过程更为复杂；此外，该微分方程的组合求解算法中耦合变量的表示方式是一个子系统的输入向量直接由另一个子系统的状态变量获得，复杂产品协同仿真的实际情况比这复杂，一般而言，某个子系统的输入向量往往来自另一个子系统的状态变量、输入向量和时间的函数。因此，需要对该组合算法加以扩展，以支持复杂产品多学科协同仿真所存在的这几种实际情况。

2 基于联合仿真步的协同仿真求解算法

对于式(1)表示的某个完整系统的协同仿真模型，分解成两个耦合的子模型式(2)、(3)，该复杂产品多学科协同仿真模型的求解过程可以用图3所示来表示。

图3 基于协调器的协同仿真模型求解

每个子模型由一组微分方程来描述，并由各自的积分器求解，每个子模型包含两组耦合向量，即输入向量和输出向量，来表示子模型之间的交互情况。协调器采用基于组合算法的思想，协调每个子模型积分器的求解过程。当推进一个积分步时每个子模型积分器向协调器发送输出向量，同时接收输入向量。为了使积分步之间能够衔接起来，每个子模型

的输入向量都采用插值算法获得，即图3中的2112,Y

Y ，其中12,I I 为插值公式。根据组合算法的运行结果，整个积分过程能获得和完整模型积分相同的效果。实际上，图3所示的模型采用的组合算法是一种基于微分方程的组合算法，它和数值计算领域求解大型常微分方程组时采用的组合算法是比较相关的。

定理1要求子系统间的数据交换在积分步上进行。然而，对于复杂产品的多学科协同仿真，其仿真步长的联合推进是由学科软件的外部步长控制，学科软件在提供外部接口时，只允许用户从一个时间段推进到另一个时间段，而内部的积分过程不可控，因此需要将定理1的条件放宽，即：使子系统间的数据交换在联合仿真步上进行，联合仿真步长一般包含多个积分步，这样就可以研究得到基于联合仿真步的组合算法。

为简化起见，我们在这里只考虑联合仿真步是等步长的

情况，来分析式(1)的分解模型式(2)、(3)的求解算法。

如图4所示，用点0,0,1,2,,i t t iH i N =+="将0[,]N t t 离散化成N 等份，H 为联合仿真步长。用点,111,0,1,2,,1i j ij j i t t h j M +=+=?"将区间1[,]i i t t +离散化成

1i M 份，0,11,i i i M i i t t t t +==，1j h 为子系统1在区间1[,]i i t t +上的内部积分步长，10111,11M i H h h h ?=+++"。类似地，用点,122,0,1,2,,1i j ij j i t t h j M +=+=?"将区间1[,]i i t t +离散化成

2i M 份，0,21,i i i M i i t t t t +==，2j h 为子系统2在区间1[,]i i t t +上的内部积分步长，20212,21M i H h h h ?=+++"。

图4 基于联合仿真步的等步长推进过程

设已得到1221,Y Y 的计算值序列1221,,0,1,,k k Y Y k i ="。记

01(,,,)T i i T t t t ="，112012112(,,,)T T T T i i Y Y Y Y ="，221021121(,,,)T T T T

i i Y Y Y Y ="，

构造分别由1(,)i i Y T 和2(,)i i Y T 确定的Lagrange 插值函数11(,)()Y i i i I Y T t 和22(,)()Y i i i I Y T t 。

令2122()(,)()i Y i i i Y

t I Y T t = ，代入(2)得到初值问题： 1112111()((),(),),()i i i dz t f z t Y

t t z t z dt

== (7) 用单步积分公式以步长1j h 数值积分，得到递推式：

1,1111121111011,111

(,,,)0,1,,1,,i j ij j f ij i ij j i i i i M i i z z h z Y

t h j M z z z z ψ++=+=?== " (8)

其中，增量函数1121(,,,)f z Y t h ψ是由所采用的单步积分公式

和1f 确定的，仅是1211,,,ij i ij j z Y

t h 的函数。求出11i z +后，按照输出公式121121()((),(),)Y t y z t Y t t =代入

11i z +，211()i i Y

t + 和1i t +求出121()i Y t +即121i Y +。令12111111()(,)()i Y i i i Y

t I Y T t ++++= ，代入(3)得到初值问题： 22212122()((),(),),()i i i dz t f z t Y t t z t z dt

+== (9) 用单步积分公式以步长2j h 数值积分，得到递推式：

2,12222121222022,221

(,,,)0,1,,1,,i j ij j f ij i ij j i i i i M i i z z h z Y

t h j M z z z z ψ+++=+=?== " (10)

其中，增量函数2212(,,,)f z Y t h ψ是由所采用的单步积分公式

和2f 确定的，仅是21212,,,ij i ij j z Y

t h + 的函数。求出21i z +后，按照输出公式212212()((),(),)Y t y z t Y t t =代入

21i z +，121i Y +和1i t +求出211()i Y t +即211i Y +。

这样就描述了一个从i t 推进到1i t +的完整过程。将公式

(8)(10)中的i 用1i +代替，不断重复这个过程，直到i N t t =为止，计算过程如图5所示。我们称这样的算法模型为基于联合仿真接口的等步长组合算法。在每个联合仿真步上，计算的误差将比在每个积分步上更大，因此要使该组合算法收敛，联合仿真步必须取得小一些。

协调器

组合算法

Y 12

Y 21

Y 12

Y 21

子模型积分器

()11121,,z

f z Y t = ()1

1211212121,,,I Y y z Y

t Y Y =→ ()22212,,z

f z Y t = ()2

2122121212,,,I Y y z Y t Y Y =

→ h 1

h 2

t i

t i+1

S 1S 2

2009年10月梁思率，

等：复杂产品多学科模型的协同仿真求解算法研究 Oct., 2009

图5 等步长组合算法串行计算流程图

在上述算法模型中，一个完整步的计算次序为：先构造21()i Y

t (特别地，当插值阶为0时2121()i i Y t Y = )，进行若干次外插和数值积分由(7)式得到11i z +，代入输出公式(2)求出

121i Y +；再构造121()i Y

t + (特别地，当插值阶为0时121121()i i Y t Y ++= )，进行若干次内插和数值积分(9)式得到21i z +，代入输出公式(3)求出211i Y +。该求解算法中所涉及的插值方法和数值积分方法，可以根据实际需要进行灵活选取。

3 算法应用与实例

在产品设计领域，ADAMS 是应用非常广泛的动力学分析软件，MATLAB 常用于控制领域的建模计算，目前

ADAMS 提供了与MATLAB 的商用软件接口，支持集中式的协同仿真。我们通过对该软件的分析认为，商用仿真接口提供的离散求解模式是一种基于联合仿真步的组合算法。以

ADAMS/CONTROL 和MATLAB/SIMULINK 的接口为例，其协同交互如图6所示，让我们来具体描述利用商用联合仿真接口求解多学科系统模型的过程：

(1) 首先，

用ADAMS 求解器对ADAMS 模型进行计算，从初始的i t 时刻推进一个仿真步长H 至1i i t t H +=+时刻，利用MATLAB 模型i t 时刻的输出作为ADAMS 模型i t 时刻的输入。这时ADAMS 求解器采用内部的积分方法将模型的状态从i t 时刻开始积分，经过若干个积分步得到1i t +时刻的状态，并计算1i t +时刻的输出。

(2) 其次，

用MATLAB 求解器对MATLAB 模型进行计算，从初始的i t 时刻推进一个仿真步长H 至1i i t t H +=+时刻，利用ADAMS 模型1i t +时刻的输出作为MATLAB 模型i t 时刻的输入。这时MATLAB 求解器采用内部的积分方法将模型的状态从i t 时刻开始积分，经过若干个积分步得到1

i t +时刻的状态，并计算1i t +时刻的输出。

(3)用1i +代替i ，重复(1)、(2)过程，直到仿真时间推进

到终止时间为止。

图6 基于联合仿真接口的数据交互过程

图7 基于联合仿真接口的组合算法

容易看出，它的插值阶为0，且子系统间的速率比为1即推进步长相同(当然，子系统间的推进步长也可以呈整数倍的关系)。因此，该算法是等速率等步长推进算法。这种算法由于插值是采用包含最新值的插值多项式而成为串行算法。

为验证这种组合算法的有效性，我们选取某航天产品模型作为多学科协同仿真实例。该模型由动力学、控制、气动三个子模型组成。分别建立每个子模型对应的联邦成员，在

HLA 环境下运行仿真。设仿真总时间为2秒，联合仿真步长为0.005秒，

得到仿真结果如图8所示(因仿真对象的限制，图中省略了所表示的具体仿真参数)。

4030201000.5

1.5

21.6

1.41.20.810.60.40.20.5 1 1.5

2t i t i

t i-1t i+1t i+1

t i-1

ADAMS

MATLAB

2009年10月系统仿真学报 Oct., 2009

图8 多学科协同仿真结果

4 结论

复杂产品的多学科协同仿真是面向一个多领域模型的耦合系统，通过多个属于不同学科的子模型之间的并行运行、实时交互，共同实现对涉及学科间交互的复杂产品设计问题的仿真分析，其仿真原理都是将每个子模型交由各自的仿真软件进行解算，子模型之间则采用一定的交互机制完成输入输出数据的更新。产品开发领域的多学科协同仿真通常是连续系统仿真，其本质是微分方程的数值积分过程。类似于数值计算中采用分解算法求解常微分方程组，多学科协同仿真也是采用基于多个求解器的组合算法求解产品模型。本文采用数值计算领域求解大型常微分方程组的分解算法，研究了多学科协同仿真的求解算法，并通过一个具体实例验证了该算法的有效性和收敛性。

参考文献:

[1]

王克明. 多学科协同仿真平台研究及其应用[D]. 清华大学, 2005.

[2]

沈俊, 宋健. 基于ADAMS 和Simulink 联合仿真的ABS 控制算法研究[J]. 系统仿真学报, 2007, 19(5): 1141-1143. (Shen Jun, Song

Jian. ADAMS and Simulink Co-Simulation for ABS Control Algorithm Research [J]. Journal of System Simulation (S1004-731X), 2007, 19(5): 1141-1143.)

[3]

黄先祥, 马长林, 高钦和, 等. 大型装置起竖系统协同仿真研究

[J]. 系统仿真学报, 2007, 19(1): 1-2. (Huang Xian-xiang, Ma Chang-lin, Gao Qin-he, et al. Studies for Collaborative Simulation of Large Erection Mechanism System [J]. Journal of System Simulation (S1004-731X), 2007, 19(1): 1-2.) [4]

陈晓波, 熊光楞, 郭斌, 等. 基于HLA 的协同仿真运行研究[J]. 系统仿真学报, 2003, 15(12): 1707-1711. (Chen Xiao-bo, Xiong

Guang-leng, Guo Bin, et al. Research on Co-simulation Running Based on HLA [J]. Journal of System Simulation (S1004-731X), 2003, 15(12): 1707-1711.)

[5]

郭斌, 范文慧, 熊光楞. 基于HLA 的复杂产品多联邦协同仿真运行管理研究[J]. 系统仿真学报, 2006, 18(5): 1212-1216. (Guo Bin,

Fan Wen-hui, Xiong Guang-leng. Research of HLA-Based Multi-Federation Collaborative Simulation Management System for Complex Product Design [J]. Journal of System Simulation (S1004-731X), 2006, 18(5): 1212-1216.)

[6] 刘德贵, 费景高. 动力学系统数字仿真算法[M]. 北京: 科学出版社, 2000.

[7] 费景高. 应用插值的数值求解常微分方程组初值问题的分解算法的收敛性和收敛阶[J]. 数值计算与计算机应用, 1984, (12):

209-218. [8]

宋健, 张越今. ADAMS 软件应用中解决数值发散的技巧[J]. 设计计算研究, 1996, (12): 7-11.

(上接第6367页)

5 结论

(1) 从图3看出，测头支杆相对于动平台坐标系{O p }XY

平面的垂直度误差hd zx ε、hd

yz ε对支杆伸缩误差影响极大，另外，第5、6、13、18项误差Td iyz ε、()x i εψ、()y i ε?、Tu izx ε对

支杆伸缩误差影响也较突出，这表明：

保证测头支杆与动平台有很好的垂直度，再者保证虎克铰的制造和装备精度，可以有效提高并联坐标测量机的工件测量精度。

(2) 从图4、5看出，从支杆的伸缩误差(传感器检测环节误差)到工件测量误差不是简单的累加关系，总体上有一定的弱化，尤其在与基础平台平行的方向上测量误差更加减弱，这表明：利用并联坐标测量机测量工件时，基础平台平面不要正对着工件，即要使得工件与基础平台中心连线与基础平台平面成锐角，以提高测量精度。

(3) 本文运用所提出的构件、副误差矩阵分析了6-TPS 型并联坐标测量机的误差因素，建立了工件测量误差模型，也为误差辨识与补偿奠定了基础。

参考文献:

[1]

LIU Dejun, CHE Rensheng, HUANG Qingcheng, et al . An error-

model study of a new coordinate-measuring machine based on the parallel-link mechanism [J]. Measurement Science and Technology (S0957-0233), 2000, 11(12): 1721-1725.

[2]

WANG S M, EHMANN K F. Error model and accuracy analysis of a six-DOF Stewart platform [C]// Proc. of the 1995 ASME International Mechanical Engineering Congress and Exposition, San Francisco, CA, USA, 1995, 2(1): 519-530.

[3]

孟婥, 车仁生. 并联六坐标测量机的误差模型和误差补偿[J]. 哈尔滨工业大学学报, 2004, 36(3): 317-320.

[4] 陈修龙, 赵永生, 邓昱. 新型并联机器人坐标测量机误差建模与仿真[J]. 计算机集成制造系统, 2008, 14(3): 477-481.

[5] 朱小蓉, 沈惠平. 一种新颖并联运动坐标测量机的误差建模与仿真[J]. 中国制造业信息化, 2006, 35(17): 56-60.

[6] 张秀峰, 季林红, 孙立宁. 精密并联机器人系统误差的分析与补偿[J]. 清华大学学报, 2005, 45(5): 630-633.

[7] 刘又午, 章青, 赵小松, 等. 数控机床全误差模型和误差补偿技术的研究[J]. 制造技术与机床, 2003, (7): 46-50.

[8] 刘又午, 章青, 赵小松, 等. 基于多体理论模型的加工中心热误差补偿技术[J]. 机械工程学报, 2002, 38(1): 127-130.

[9] 刘丽冰, 王广彦, 刘又午. 复杂机械系统运动误差自动建模技术研究[J]. 中国机械工程, 2000, 11(6): 642-646.

[10] 刘延斌, 韩秀英. 一种五轴数控机床的综合误差建模与补偿[J].

制造技术与机床, 2007, (11): 35-38.

2×104

1.50.5100.5

1.5

2 -4

-3-2-1010.5

1.5