【技术干货】GCkontrol连续求解器在建模仿真中的应用
本文介绍了GCKontrol提供的离散求解器和连续求解器的基本理论和使用方法，通过示例展示了离散求解器和连续求解器在建模仿真中的应用。
微分方程：表示未知函数和未知函数的导数以及自变量之间的关系的方程，叫做微分方程。
常微分方程（Ordinary Differential Equations，ODE）：如果在一个微分方程中出现的未知函数只含一个自变量，这个方程叫做常微分方程。
定义：如果对于
，只要
满足ODE方程且
，就有
对所有的
成立，则称ODE方程
的解是稳定的。
下面用三个例子说明ODE解的稳定性。
渐进稳定解：
图1 渐进稳定
稳定但不是渐进稳定：
图2 稳定但不是渐进稳定
不稳定：
图3 不稳定
总结：ODE方程
的解为：
。如果
：指数发散：每一个解都是不稳定的；如果
：每一个解都是渐进稳定的。
一些关键核心技术实现破，战略性新兴产业发展壮大，载人航天、探月探火、深海深地探测、超级计算机、卫星导航、量子信息、核电技术、大飞机制造、生物医药等取得重大成果，进入创新型国家行列。
通常要仿真动态系统，可以通过计算其在指定时间跨度内连续时间步的状态来实现。时间步是发生计算的时间间隔，时间步的大小称为步长。以这种方式计算模型状态的过程称为解算模型。
对于表示动态系统连续状态的ODE方程，往往通过离散化，采用数值微分法、数值积分法和泰勒逼近法的思想，建立数值解。
数学家们开发出了多种数值积分方法来解算ODE方程，包括欧拉（Euler）方法、龙格-库塔（Runge–Kutta）方法等。基于这些数值积分方法，GCKontrol提供了一组程序用来解算模型，称为求解器。每个求解器代表一种特定的模型解算方法。
GCKontrol中的求解器分为定步长离散求解器和定步长连续求解器。定步长求解器使用固定的步长从仿真开始到仿真结束来解算模型。
当使用离散求解器进行求解时，首先将模型离散化（生成对应的离散代码），然后通过欧拉方法迭代计算模型的下一个仿真时间步的状态。
当使用连续求解器进行求解时，模型和求解器是分开的。首先对模型计算状态变化率（dU/dt），求解器积分这些dU/dt，根据模型在上一个时间步的状态和状态导数来计算模型在当前时间步的连续状态。连续求解器可以有不同的求解方法，在GCKontrol中采用的是四阶的龙格-库塔方法。
在编译和仿真模型时，可以基于以下因素选择合适的求解器：
系统动态特性
嵌入式代码生成
模拟离散/连续（物理）系统特性
求解器稳健性
GCKontrol支持嵌入式代码生成，GCKontrol工程在编译代码时，会将工程中的模块信息、参数信息以及状态图信息等生成满足合规性标准的C/C++代码。生成的代码可以在其他C/C++集成开发环境下运行，进行代码编译和调试功能，用于后续的嵌入式仿真需求。
因为架构和函数不一样，离散求解器和连续求解器生成的代码有区别。实际生产生活中，绝大部分控制器的工作是离散的，对于控制算法开发来说，离散求解器更适合做嵌入式代码生成。所以在用GCKontrol做控制系统建模仿真时，被控对象模型用连续求解器求解，控制算法模型用离散求解器求解是比较通用的做法。
GCKontrol中的离散求解器是用一阶精度的欧拉（Euler）方法求解的，当选择连续求解器时，用四阶的龙格-库塔（Runge–Kutta）方法求解。
ODE方程
的解为
欧拉迭代公式
即
是增长因子，如果
，欧拉方法是稳定的，如果
，欧拉方法是发散的，其中，h为步长。
令
则不等式
可表示为
求解该不等式，可得z的取值范围如下图绿色 区域所示。绿色 区域即是欧拉方法求解的稳定区域。在绿色 区域内，数值解是稳定的，超出了绿色 区域，数值解是不稳定的。
图4 欧拉方法求解的稳定性
分析可得，求解器的时间步长h和系统动力学之间存在关系，系统的稳定区域和步长有关，步长越小，系统稳定区域越大。
在Re坐标轴上，z的最小值为-2，
即
如果取仿真步长为1ms，
则
对于使用欧拉方法求解的系统来说，仿真步长取1ms，则系统动力学特性高于300Hz的模型可能会发散。在实际应用中，系统动态特性为50Hz或100Hz的模型，采用1ms求解步长是合适的。
图5 欧拉方法求解的稳定性分析
对欧拉方法求解的稳定区域进一步分析可以发现，欧拉方法对于低阻尼系统的求解效果非常不好，往往难以正确求解或者导致结果发散。
连续求解器使用四阶龙格-库塔（RK4）公式计算下一时间步的模型状态，作为状态当前值和状态导数的显式函数。
龙格-库塔法的思想：利用f(x,y)在某些特殊点上的函数值的线性组合来构造高阶单步法的平均斜率。
四阶龙格-库塔公式：
其中
推导得
其中
是增长因子，如果
，龙格-库塔方法是稳定的，如果
，龙格-库塔方法是发散的。
求解该不等式，可得z的取值范围如下图黑色 区域所示。黑色 区域即是龙格-库塔方法求解的稳定区域。在黑色 区域内，数值解是稳定的，超出了黑色 区域，数值解是不稳定的。
图6 龙格-库塔方法求解的稳定性
下图是两种求解器的稳定区域对比。绿色 区域是欧拉方法求解的稳定区域，黑色 区域是龙格-库塔方法求解的稳定区域。
可以看出，龙格-库塔方法的求解稳定区域比欧拉方法的更大。对于低阻尼系统，采用龙格-库塔方法的连续求解器的适用性更强，可以更大程度上保证解的稳定和精度。
图7 两种求解器的稳定区域
在时间设置和子系统的属性设置中可选择求解器类型。
图8 时间配置中求解器的设置
时间配置的求解器选项可设置离散求解器和连续求解器。当时间配置中选定求解器后，除单独设置求解器为离散/连续求解器的子系统外，全部按照时间配置的求解器解算。
图9 子系统属性设置中求解器的设置
子系统的求解器设置分为离散求解器、连续求解器和沿用父级，默认选择沿用父级。当子系统选择离散/连续求解器时，子系统中的模块按照连续求解的方式求解，当子系统选择沿用父级后，子系统中模块的求解方式按照父级子系统的求解方式求解。
下图展示在时间配置和子系统设置中选择不同求解器类型的效果。
图10 求解器的设置
图中子系统Subsystem1和Integrator模块采用离散求解器求解，子系统Subsystem2采用连续求解器求解。
质量弹簧阻尼器
以质量弹簧阻尼器示例工程为例，分别用离散求解器和连续求解器求解位移的数值解，然后和位移的精确解进行对比，结果见下图。（详细的建模过程请查看示例工程《质量弹簧阻尼器》）
图11 求解结果对比
从图中可以看到，采用连续求解器求解的结果与精确解一致，采用离散求解器求解的结果与精确解存在偏差，说明在质量弹簧阻尼系统的仿真中，采用连续求解器的精度更高。
液压系统模型
以液压系统模型示例工程为例，采用离散求解器，修改时间配置的调度周期为0.0001，工程发散。（详细的建模过程请查看示例工程《液压系统模型》）
图12 结果发散
将工程设置为连续求解器求解，结果收敛。
图13 结果收敛
从该示例可以看出，即使ODE稳定，求解器可能是不稳定的，这和仿真步长h相关。当选择了较小的步长时，采用离散求解器可能导致结果发散，这种情况下需要采用连续求解器保证结果的收敛。
涡扇发动机
以涡扇发动机示例工程为例，涡扇发动机模型是一个复杂模型，有很多连续状态，采用连续求解器求解。控制模型则采用离散求解器求解，便于生成合适的嵌入式代码。（详细的建模过程请查看示例工程《涡扇发动机》）
图14 涡扇发动机模型
将离散求解器的调度周期设置为10ms，如图所示。
图15 控制器模型步长设置为10ms
发动机模型可以采用较小的步长，保证仿真结果的精度，这里将连续求解器的步长设置为0.1ms。
图16 发动机模型步长设置为0.1ms
仿真时长为100ms，对比采用离散求解器的燃油流量和采用连续求解器的总推力，结果见下图。
图17 结果对比
可以看到，燃油流量的变化是离散的，总推力的变化是连续的。这是采用不同的求解器所导致的差异，这种差异是符合实际情况的。
GCKontrol支持对模型仿真和代码生成期间模型状态的求解设置离散求解器和连续求解器。离散求解方法稳定区域较窄，对大型复杂模型进行建模时，需要保证解的数值稳定性，此时需选择连续求解器来提高求解器的阶数，增加解的稳定性区域。对于控制系统建模仿真，离散求解器更适合做嵌入式代码生成，通常将离散求解器用于控制算法模型，被控对象则采用连续求解器，提高解的稳定性和精度。
世冠科技成立于2003年，是一家专业从事工业软件系统仿真技术开发与应用的国家级高新技术企业，北京市企业科技研究开发机构、专精特新企业。为复杂装备研制单位和工业制造企业，提供可支撑产品设计研发及使用运维、覆盖产品全生命周期的完全自主研发的系统仿真工业软件和数字孪生解决方案，是国产系统仿真软件领域的领军者。
