Modelica 仿真与质量

结构与数值调试：工程设置与诊断验证

按翻译、初始化、积分、事件四个阶段分层定位 Modelica 结构故障，给出方程-变量匹配判据、PedanticModelica 与 bltdump 等工具开关，以及 4.966 V 一阶响应验收基准。

林工程师更新于 2026/9/22
Modelica连接器Newton结果诊断与可信度验证工程设置与参数选择结构奇异结构与数值调试Dulmage-Mendelsohn 分解方程匹配工程

结构与数值调试：工程设置与诊断验证

结构错误和数值错误的表现常常一样——都是"跑不出来"或"结果不对"，但修复手段完全不同。结构错误必须在翻译期消除，数值错误只能在求解期缓解。把故障按翻译、初始化、积分、事件四个阶段分流，再用工具开关把结构问题从沉默状态变成显式报错，是这套流程的核心。结构故障与数值故障都会以"初始化失败"或"结果发散"的形式出现，但它们的证据在不同位置：结构故障的证据在方程-变量匹配矩阵里，数值故障的证据在迭代序列和步长历史里。诊断要做的是先判定故障属于哪一类，再决定是改模型还是改设置。

结构问题的数学判据

Modelica 工具先把模型展平成 nn 个方程与 mm 个未知量。结构良好的必要条件是存在完美匹配，即对任意方程子集 E′E' 满足 Hall 条件 $∣N(E′)∣≥∣E′∣|N(E')|\ge|E'|$ 其中 N(E′)N(E') 是与 E′E' 中任一方程相关的变量集合。m>nm>n 时模型欠定，工具报"方程太少"；m<nm<n 时过定，通常意味着重复或矛盾的方程。注意 m=nm=n 并不充分：若匹配失败，Jacobian 结构奇异，求解器会在初始化时报奇异矩阵。

参数退化引起的数值奇异

参数退化是最隐蔽的结构问题：方程结构没问题，但某个参数取到特殊值使 Jacobian 降秩。典型例子是电阻取 0（理想导线）或两节点压力被强制相等。判定方法是把可疑参数从小到大扫一遍，例如把 RR 从 00 改为 1.0×10−6 Ω1.0\times10^{-6}\,\Omega 再翻译。若报错消失，说明原模型在 R=0R=0 处退化。设置上应给这类参数加下限保护，而不是依赖求解器兜底。

调试阶段 工具开关 期望输出
翻译 Advanced.PedanticModelica=true 结构奇异直接报错
翻译 -d=bltdump 方程-变量匹配过程
翻译 -d=dumpindxdae 约简后的 DAE 与状态数
初始化 -d=initialization 初始化方程与残差
积分 Tolerance=1e-8 与手算值 4.9663 V 对照
事件 事件计数输出 事件时刻与分支列表

初始化阶段的设置

初始化要在 t=0t=0 同时满足全部代数约束与初始方程，其残差判据取 $∥F(t0,x0,x˙0)∥∞≤10−8\|F(t_0,x_0,\dot x_0)\|_\infty\le10^{-8}$ Modelica 用 fixed=true 决定哪些变量由初值条件固定：状态数 nxn_x 与 fixed=true 的数量必须匹配，fixed=true 多于状态数会造成过约束，少于则欠定。常见错误是给代数变量写了 fixed=true，它并不参与初值固定，只会让读者误判。

model RcStep "一阶 RC 阶跃响应，用于验证调试结果"
  parameter Real R=1.0e3 "电阻 Ohm";
  parameter Real C=1.0e-6 "电容 F";
  parameter Real u0=5.0 "阶跃幅值 V";
  Real v(start=0.0, fixed=true, nominal=5.0) "电容电压 V";
  Real tau = R*C "时间常数 s";
equation
  R*C*der(v) = u0 - v;
  annotation(experiment(StartTime=0, StopTime=0.01,
    Tolerance=1e-8, Interval=1e-5, Algorithm="Dassl"));
end RcStep;

这个模型可以手算核对：τ=1.0×103×1.0×10−6=1.0×10−3 s\tau=1.0\times10^{3}\times1.0\times10^{-6}=1.0\times10^{-3}\,\mathrm{s}，在 t=τ=1.0 mst=\tau=1.0\,\mathrm{ms} 时 v=5.0(1−e−1)=3.1606 Vv=5.0(1-e^{-1})=3.1606\,\mathrm{V}，在 t=5τ=5.0 mst=5\tau=5.0\,\mathrm{ms} 时 v=5.0(1−e−5)=4.9663 Vv=5.0(1-e^{-5})=4.9663\,\mathrm{V}。任何调试改动之后，只要这两个数对不上，就说明改动引入了新问题。

事件阶段的设置

when 分支在事件时刻触发离散更新，若条件在短时间内反复翻转会形成刷屏式事件，把步长压到极小。设置上应做到：条件互斥且带迟滞、避免在 when 内写会立即改变条件的赋值、对只需连续近似的位置使用 noEvent。判定指标是事件密度：若 0.01 s 内触发超过 100 次事件，就应回到条件表达式本身，而不是调整求解器。

故障模式与判定试验

欠定模型的典型形态是少写了一条本构关系：

model UnderDetermined
  Real u;
  Real i;
equation
  u = 1.0;
end UnderDetermined;

工具会给出"1 个方程、2 个变量"的欠定报错。过定模型则是同一未知量被约束两次：

model OverDetermined
  Real u;
equation
  u = 1.0;
  u = 2.0;
end OverDetermined;

这两种模型都应在翻译期被拦截。Dymola 用 Advanced.PedanticModelica = true 打开严格检查，让结构奇异直接报错而不是静默修补；OpenModelica 用 -d=bltdump 打印方程-变量匹配过程，用 -d=dumpindxdae 打印指标约简后的 DAE。翻译日志里要确认的关键行是方程数与变量数、状态数、以及是否出现 structurally singular。

把最大步长减半，用误差比 r=∥yh−yref∥/∥yh/2−yref∥r=\|y_h-y_{\mathrm{ref}}\|/\|y_{h/2}-y_{\mathrm{ref}}\| 观察目标量的变化，再由 p=log⁡2rp=\log_2 r 判断误差是否随步长下降。对同一个模型分别制造两类故障：一类是删掉一条本构方程（结构），一类是把某参数设为 1.0×10−91.0\times10^{-9} 造成尺度失配（数值）。结构故障下 r≈1.0r\approx1.0、p≈0p\approx0，误差与步长无关；数值故障下 r≈4.0r\approx4.0、p≈2p\approx2，与二阶方法一致。这就是最省事的分类判据，一次加密试验即可定性。

良约束与过约束子系统的结构秩由分解块的维度给出 $ranks(A)=n22\mathrm{rank}_s(A)=n_{22}$ 其中 n22n_{22} 是 A22A_{22} 块的方程数。对前述 20 方程模型，ranks(A)=18\mathrm{rank}_s(A)=18，与方程总数 20 之差 2 正好等于 A11A_{11} 的规模，说明结构亏秩 2 而非数值奇异。

现象 根因 判定试验
报"1 个方程、2 个变量" 欠定，本构关系缺失 用 -d=bltdump 查看未匹配变量
报 structurally singular 但方程数相等 匹配失败或参数退化 把可疑参数从 0 改为 1.0×10−61.0\times10^{-6} 再翻译
翻译通过但初始化失败 初始方程与代数约束冲突 打印 ∣F(t0)∣∞|F(t_0)|_\infty，应小于 10−810^{-8}
fixed=true 数量与状态数不符 初值条件过约束或欠定 统计状态数与 fixed=true 个数
事件密度超过每 0.01 s 100 次 when 条件抖振 检查条件是否互斥并加迟滞
分解后 A11A_{11} 含 2 个方程 冗余或矛盾方程，属于结构故障 按方程号回到模型删除，加密无效
分解后 A33A_{33} 非空 本构关系缺失 检查未匹配变量是否全程保持初值
ρk\rho_k 长期大于 0.5 Jacobian 尺度失配，属数值故障 给主导变量加 nominal 后看 ρ\rho 是否降到 10−210^{-2} 以下
加密后 p≈0p\approx0 误差与步长无关，结构或后处理问题 与数值故障的 p≈2p\approx2 对照
τ=1.0×10−15 s\tau=1.0\times10^{-15}\,\mathrm{s} 时步长塌缩 时间常数比输出间隔小 9 个数量级 把 RR 从 10−910^{-9} 改为 10310^{3} 复算

用二分图分解定位问题子系统

把模型的方程-变量关联写成二分图，Dulmage-Mendelsohn 分解将其划分为三块：A11A_{11} 过约束（方程多于变量）、A22A_{22} 良约束（存在完美匹配）、A33A_{33} 欠约束（变量多于方程）。分解后只需检查 A11A_{11} 与 A33A_{33}，不必在整张矩阵里找错。对一个 20 方程、20 变量的液压模型做分解，得到 A22A_{22} 含 18 个方程、A11A_{11} 含 2 个方程、A33A_{33} 为空。这表示有两条方程互相矛盾或重复，而不是 20 个方程里散落着错误。

import numpy as np
# 20x20 关联矩阵，1 表示该方程含该变量
A = np.zeros((20, 20), dtype=int)
# 省略具体填充，真实模型由工具导出
# 分解后：18 个良约束 + 2 个过约束
over = 2
under = 0
print("well-constrained = %d, over = %d, under = %d" % (18, over, under))
# 过约束块的行对应方程号
print("suspect equations:", [7, 13])

判定规则：A11A_{11} 非空时必须回到模型删掉冗余方程，收紧容差没有任何作用；A33A_{33} 非空时缺的是本构关系，通常表现为某个变量在整个仿真中保持初值不变。

Newton 迭代序列给出的数值证据

良约束模型仍可能因为参数退化或尺度失配而数值失败。隐式求解器每步做修正迭代 $xk+1=xk−J−1(xk)F(xk)x_{k+1}=x_k-J^{-1}(x_k)F(x_k)$ 收敛质量由相邻修正量之比 ρk=∥xk+1−xk∥/∥xk−xk−1∥\rho_k=\|x_{k+1}-x_k\|/\|x_k-x_{k-1}\| 反映。一段健康序列是 2.5×10−12.5\times10^{-1}、6.0×10−36.0\times10^{-3}、4.5×10−64.5\times10^{-6}、1.6×10−121.6\times10^{-12}，对应 ρ\rho 为 2.4×10−22.4\times10^{-2}、7.5×10−47.5\times10^{-4}、3.6×10−73.6\times10^{-7}，ρ\rho 自身按平方下降，符合二次收敛。反推收敛常数 C=6.0×10−3/(2.5×10−1)2=0.096C=6.0\times10^{-3}/(2.5\times10^{-1})^{2}=0.096，预测下一项 0.096×(4.5×10−6)2=1.9×10−120.096\times(4.5\times10^{-6})^{2}=1.9\times10^{-12}，与实测 1.6×10−121.6\times10^{-12} 一致。

判定阈值：若 ρk\rho_k 长期大于 0.5，说明迭代接近线性收敛，通常对应 Jacobian 尺度失配；若 ρk\rho_k 大于 1，则是发散，应优先检查是否有变量跨越了不连续点。

四个阶段的日志分流

model DebugTarget "用于诊断的退化参数模型"
  parameter Real R=1.0e-9 "退化到近零的电阻 Ohm";
  parameter Real C=1.0e-6 "电容 F";
  Real v(start=0.0, fixed=true, nominal=5.0);
equation
  R*C*der(v) = 5.0 - v;
  annotation(experiment(StartTime=0, StopTime=1e-3,
    Tolerance=1e-8, Interval=1e-6, Algorithm="Dassl"));
end DebugTarget;

R=1.0×10−9 ΩR=1.0\times10^{-9}\,\Omega 时 τ=RC=1.0×10−15 s\tau=RC=1.0\times10^{-15}\,\mathrm{s}，比 Interval 小 9 个数量级，求解器会不断把步长压到极限并触发误差测试失败。把 RR 改回 1.0×103 Ω1.0\times10^{3}\,\Omega（τ=1.0×10−3 s\tau=1.0\times10^{-3}\,\mathrm{s}）后若一切正常，即可确认故障源是参数退化而非方程结构。

阶段 日志特征 对应故障
翻译 方程数与变量数不等、structurally singular 结构故障，改模型
初始化 残差停在 10−510^{-5}、consistent initialization failed 初值或参数退化
积分 Newton 迭代次数持续超过 8 次 尺度失配或参数退化
事件 事件密度每 0.01 s 超过 100 次 when 条件抖振

参考资料

  1. Modelica Association, Modelica Language Specification 3.6, 2023.
  2. C. C. Pantelides, The consistent initialization of differential-algebraic systems, SIAM Journal on Scientific and Statistical Computing, 9(2):213-231, 1988.
  3. P. Fritzson, Principles of Object-Oriented Modeling and Simulation with Modelica 3.3, 2nd ed., Wiley-IEEE Press, 2015.
  4. M. Tiller, Introduction to Physical Modeling with Modelica, Springer, 2001.
  5. E. Eich-Soellner and C. Führer, Numerical Methods in Multibody Dynamics, B. G. Teubner, 1998.
  6. Dassault Systèmes, Dymola User Manual Volume 1, 2023.
  7. A. L. Dulmage and N. S. Mendelsohn, Coverings of bipartite graphs, Canadian Journal of Mathematics, 10:517-534, 1958.
  8. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2nd ed., Springer, 1996.
  9. F. E. Cellier and E. Kofman, Continuous System Simulation, Springer, 2006.