Co-Simulation Master Algorithms: Jacobi vs Gauss-Seidel, SSP and Algebraic Loops

Export two models as Functional Mock-up Units (FMUs), wire one’s output to the other’s input, and run both at a 100 ms communication step. Each FMU is individually validated and each solver is individually accurate, yet the coupled result can drift, ring, or blow up. In the toy system built later in this post, the very same pair of subsystems produces a worst-case position error of 51 (against a motion amplitude of 1) under one stepping scheme and 0.24 under another, with identical step size and identical subsystem solvers. The difference is entirely the co-simulation master algorithm: the orchestration logic that decides who steps when, what inputs each unit sees, and what happens when a step turns out to be wrong.

This matters in 2026 because digital twins increasingly compose vendor models, controller code, and physics from different tools, and the composition layer is where credibility is lost. You leave with a mental model of Jacobi versus Gauss-Seidel stepping, a way to reason about coupling error, rules for step-size control with rollback, tactics for algebraic loops, and a runnable Python master with measured error numbers.

What this covers: the FMI co-simulation contract, serial and parallel schemes, coupling error and extrapolation, rollback, fixed-point and Anderson iteration for loops, SSP packaging, tooling, real-time and hardware-in-the-loop (HIL) operation, and a decision checklist.

Context and Background

Co-simulation couples independently developed simulators that each carry their own numerical solver. Instead of one integrator seeing one big equation system (the monolithic approach, as in a single Modelica compilation), each subsystem integrates itself and the units exchange data only at discrete communication points. This is the premise of the Functional Mock-up Interface (FMI), the Modelica Association standard that packages a simulator as an FMU with a defined C API. The motivation is practical: suppliers will ship an FMU but never their equations, and a plant model written in one tool must be driven by a controller written in another.

FMI defines two FMU kinds. Model Exchange FMUs expose derivatives and let the importer’s solver integrate them. Co-Simulation FMUs embed their own solver and advance on request, which is the interface a master algorithm drives. FMI 3.0 was released on 10 May 2022 per the specification page, and its co-simulation additions, as the FMI project’s own FAQ summarises, include early return and event mode, clocks, intermediate variable access, array variables, and a binary data type. The specification is at fmi-standard.org. Crucially, FMI standardises the FMU interface and deliberately does not standardise the master. The specification describes what a master may do (set inputs, call do-step, get and set state), not which algorithm is any good. That gap is why two compliant tools can give different answers for the same FMUs.

The other half of the picture is how a composed system is described. Before the System Structure and Parameterization (SSP) standard, every tool kept its own connection file, so a topology could not move between vendors. SSP, maintained by the Modelica Association, describes components, connectors, connections, and parameter sets in a portable archive. SSP 2.0 was announced in December 2024 (the modelica.org announcement is dated January 2025) with architecture exchange, support for FMI 3.0 based features and virtual ECUs, and more efficient parameter updates.

If you are choosing between a twin that replays recorded data and one that runs live coupled physics, read our digital twin versus simulation architecture decision guide first; co-simulation is the machinery behind the second kind. If you are still deciding how to author the component models themselves, the Modelica versus Simulink comparison covers that earlier stage, and the FMI 3.0 co-simulation FMPy tutorial shows the hands-on loading of FMUs that this post assumes.

My thesis, which you will not find in the introductory material: the master algorithm is a numerical method with an order, a stability region, and an error constant, and it should be selected and verified like one, not treated as plumbing. A co-simulation that holds inputs constant between communication points is first-order accurate in the coupling variables no matter how high-order the individual solvers are, and teams rarely measure it.

The Anatomy of a Master Algorithm

A co-simulation master algorithm is the loop that advances a set of Co-Simulation FMUs from one communication point to the next. At each macro step it exchanges outputs and inputs, calls the FMUs’ do-step function with a step size h, and then accepts the step, retries it, or shrinks it. Its schedule determines accuracy, stability, parallelism, and real-time feasibility.

Figure 1: Reference architecture. The SSP package defines topology and values; the master owns time, data exchange, state snapshots and the error estimator.

The diagram shows five responsibilities that every serious master implements, whether or not they are separated in code. The first is topology: which output feeds which input, loaded from an SSP or a tool-specific file. The second is scheduling: the order and concurrency of do-step calls. The third is data conditioning: unit conversion, input extrapolation, and filtering. The fourth is supervision: error estimation, step-size control, and rollback through saved FMU states. The fifth is the clock: either as-fast-as-possible or locked to wall time for HIL.

The FMI co-simulation contract the master relies on

An FMI 3.0 Co-Simulation FMU is stepped through fmi3DoStep. From the specification text I could verify, the call takes the current communication point and a communication step size that must be greater than zero, and it can report an early return together with the last successful time it reached. If I recall the full C signature correctly (verify against the header in your FMI version), it also carries a flag telling the FMU that the importer will not restore any state older than the current point, which lets the FMU skip internal bookkeeping, plus outputs signalling that event handling is needed and that the FMU requested termination.

Two capabilities in the FMU’s model description matter most for masters. The first is canGetAndSetFMUState: when true, fmi3GetFMUState copies the internal state to an opaque handle and fmi3SetFMUState restores it, and the specification states the FMU must not modify the saved state, so one snapshot can be restored repeatedly. That is the primitive behind rollback. The second is the ability to provide output derivatives, which FMI 3.0 retains for Co-Simulation and Scheduled Execution when the FMU supports them. Output derivatives let the master extrapolate an output with a higher-order polynomial rather than a constant. The 3.0 documentation also notes that the old mechanism for setting input derivatives was removed and replaced by allowing input setting in Intermediate Update Mode, a change that matters if you port an FMI 2.0 master.

Intermediate Update Mode itself is a new FMI 3.0 idea: the FMU calls back into the importer during a long step, so the importer can read intermediate outputs and, where permitted, supply intermediate inputs, and it can signal that the step should return early. Early return in turn lets an FMU stop a step when an internal event occurs, so the master can handle the event at its true time rather than at the end of the macro step. These features narrow the gap between co-simulation and monolithic simulation, but a master must be written to exploit them. Many tools still treat an FMI 3.0 FMU as if it were FMI 2.0.

Jacobi: everyone steps from the same snapshot

In the Jacobi (parallel, non-iterative) scheme, every FMU receives inputs computed from the outputs at time t, then all FMUs step to t plus h concurrently. No unit sees any other unit’s result for the current interval. Because the do-step calls are independent, they can run on separate threads or machines, and wall-clock cost per step is the maximum, not the sum, of the unit costs.

Figure 2: Jacobi steps all FMUs from time t using old outputs; Gauss-Seidel steps FMUs in sequence so later units see fresh outputs of earlier ones.

The cost of that parallelism is a built-in one-step delay on every coupling. With constant input held between communication points, the receiving FMU effectively sees a signal that lags the true one by on average h/2, and both sides see a lag at once. For an oscillatory interaction such as two masses on a spring, the delay behaves like a spurious negative damping in the loop, which is why Jacobi can go unstable at step sizes where the monolithic model is perfectly well behaved.

Gauss-Seidel: serial, with fresher data

In Gauss-Seidel (serial) stepping, the master orders the FMUs. The first steps using old inputs; its new outputs at t plus h are then passed to the second, which steps with fresher data, and so on down the chain. Only the first unit in the order suffers the full delay, and in a feed-forward chain the scheme is effectively exact in the coupling, because each unit sees a value that was properly computed for the interval.

The costs are serialisation (wall-clock cost is the sum, not the maximum) and ordering sensitivity. For a feedback loop, somebody has to go first and see stale data, so the choice of order is a modelling decision: step the fast, stiff, or low-latency-tolerant unit last, and the unit whose inputs vary slowly first. A common heuristic is to order along the physical causality (actuator, plant, sensor), and a graph-based master can pick the order by finding strongly connected components of the connection graph and topologically sorting the components. Within a component, Gauss-Seidel still lags once around the loop.

A decision framing

Neither scheme dominates. Jacobi wins when parallel hardware is available, when units are many and loosely coupled, and when real-time determinism matters more than the last digit of accuracy. Gauss-Seidel wins when the topology is a chain or a tree, when units are few, and when the coupling is tight enough that one-step staleness hurts. In my measurements below, Gauss-Seidel roughly halved the worst-case error compared with Jacobi at equal step size for a symmetric two-body loop, but only after the Jacobi step size was small enough to be stable at all.

Coupling Error, Extrapolation, Step Control and Algebraic Loops

Inside each FMU the solver controls its own error. The master adds a second, separate error: the coupling error from exchanging data only at communication points. It is the term that dominates in most real systems, and it is the only one the master can influence.

Where coupling error comes from

Consider two subsystems with states x1 and x2, and inputs u1 = y2 and u2 = y1. Over a macro step the master must supply u1(t) and u2(t) for t in the interval, but only knows values up to time t_n. If it holds them constant (zero-order hold, ZOH), the local error in each input is O(h) times the signal’s derivative, so the local state error is O(h squared) and the global error over a fixed horizon is O(h). That is the first-order accuracy of the simplest master. If the master instead extrapolates linearly from the last two points, the input error drops to O(h squared) and the global error to O(h squared), provided the system is stable enough not to amplify errors. Quadratic extrapolation gives third order in principle, but the extrapolation polynomial amplifies noise and any discontinuity, and the stability region shrinks as order grows. In practice, linear extrapolation is the sweet spot, and it requires remembering one past output per connection.

Where an FMU offers output derivatives, the master can use a Hermite-style predictor: value plus derivative times time. This is cheaper in memory than a two-point history and works at the very first step, but the FMU must actually implement the derivative (the capability flag is optional), and it must be consistent with the output.

Step-size control through error estimation

A fixed macro step wastes effort during quiet phases and fails during transients. Variable macro steps need an error estimate that does not require the monolithic solution. Two estimators are common. The first is step doubling (Richardson style): advance one step of h, then restore the saved state and advance two steps of h/2, and compare. It needs state restore on every FMU, so it is only available when canGetAndSetFMUState is true everywhere. The second compares the extrapolated input prediction with the actually computed output at the end of the step: if the predicted and realised coupling signals differ by more than a tolerance, the step was too long. This needs no rollback to estimate, only to react, and it is cheap because it reuses the data the master already holds.

Figure 3: Step acceptance loop. Snapshot all FMUs, step with extrapolated inputs, estimate coupling error, then accept and grow or roll back and shrink.

The control law should be the same kind of controller used in ordinary ODE solvers: new step = h times safety factor times (tolerance divided by error) raised to 1/(p+1), where p is the order of the coupling scheme, clamped to a growth factor limit (say 2) and shrink limit (say 0.2). Those numbers are conventional in ODE codes rather than FMI requirements, so tune them. Two additional rules save pain: never let the macro step exceed any FMU’s declared maximum step, and snap the step to land exactly on the next scheduled external event such as a controller sample time or a setpoint change.

Rollback: what is required and what it costs

Rollback means restoring all FMUs to the state at the start of the step and re-executing with a smaller h or better inputs. In FMI it is the fmi3GetFMUState / fmi3SetFMUState pair. Three caveats are worth knowing. First, not every FMU supports it; many commercial exports set the capability to false because the tool’s internal state includes opaque handles. Second, a state snapshot has a memory and time cost proportional to the FMU’s internal state, which for a detailed thermal-fluid model can be megabytes, so snapshotting every step may dominate runtime. The master should snapshot lazily, only when the error controller is active. Third, rollback only restores the FMU, not the world: any side effect (a file written, a message published to a broker) outside the FMU repeats. A twin that publishes results to MQTT must buffer them until a step is accepted.

The 3.0 documentation describes the flag that tells an FMU the importer will never restore a state before the current point; this is the master’s promise that lets the FMU free history. A master that does variable-step with rollback must not set it, and one that does fixed-step without rollback should.

Not all engines implement rollback. According to the INTO-CPS Maestro 2 documentation, the fixed-step Jacobi algorithm there supports rollback, while its variable-step algorithm is listed without rollback, and the documentation covers FMI 2.0 co-simulation FMUs. That is a good example of why you check what your tool actually does.

Algebraic loops

An algebraic loop occurs when an FMU’s output depends on its input within the same instant (direct feedthrough) and the connection graph closes a cycle through such dependencies. Unit A computes y_A = f(u_A), unit B computes y_B = g(u_B), and u_A = y_B, u_B = y_A. The values at a single time point must satisfy a coupled nonlinear equation, but neither FMU can be asked to solve it for you, because FMI co-simulation only lets you set inputs and read outputs, not solve for them.

FMI 3.0 FMUs declare which outputs depend on which inputs in the model description (the structure information, with dependency kinds such as dependent and constant). A good master reads this before running and builds the dependency graph. Only the cycles with direct feedthrough on every edge are true algebraic loops. A cycle that passes through at least one state (integrator-like delay) is not an algebraic loop but an ordinary feedback and is handled by stepping.

For genuine loops there are three standard remedies. The first is to break the loop by inserting an explicit delay, a one-step lag, and accept the error. The second is to iterate: guess the loop variables, evaluate the FMUs, compare, update the guess, repeat until the residual is small. Since FMI lets you set inputs and get outputs at a fixed time without stepping, the master can iterate at the initialization stage and at event instants. Within a step, iterating requires rollback, because a do-step has already advanced time. The third is to restructure the model so one side exposes the implicit equation, which is the best long-term fix and often impossible with supplier black boxes.

Plain fixed-point (Picard) iteration u_new = G(u) converges when the loop map is a contraction, roughly when the product of loop gains has magnitude below 1. When it is above 1, the iteration diverges, and when it is close to 1 it is very slow. Anderson acceleration keeps the last m iterates and residuals and chooses the next guess as the combination that minimises the residual in a least-squares sense. For a scalar affine loop, Anderson with m = 1 is equivalent to the secant method and solves it exactly in a few evaluations.

I checked this with a small script (included in the examples below). The two-unit loop is y1 = a u1 + 1, y2 = b u2 – 2, with u1 = y2 and u2 = y1, and tolerance 1e-10 on the iterate change. I measured:

Loop gain a*b	Plain fixed-point	Anderson m=1
0.5	converged in 35 iterations	3 iterations
0.9	did not converge in 200 iterations	3 iterations
-1.5	diverged	3 iterations
2.0	diverged	3 iterations

The affine setting is the friendliest possible, so treat the Anderson numbers as a lower bound on difficulty; real loops are nonlinear, and Anderson then needs globalisation (damping, safeguards) and typically takes more iterations. The qualitative lesson holds: the convergence criterion of plain iteration is a property of the model, not the master, and Anderson or Newton-type updates remove that dependence at the cost of more state in the master.

A Runnable Toy Master: Two Coupled Masses

Theory is cheap, so here is a complete, runnable master that compares the schemes on a system where the monolithic exact solution is known. Mass 1 (1 kg) and mass 2 (2 kg) each hang on their own wall spring (20 and 30 N/m) and are connected to each other by a spring of 40 N/m with a small damper of 0.2 N·s/m. Mass 1 starts displaced by 1 m; everything else starts at rest. The system is linear, so the exact monolithic response comes from a matrix exponential of the full 4-state system. Each “FMU” is one mass; it integrates its own two states exactly for a given input signal (via an augmented matrix exponential), which means the only error in the experiment is coupling error. That isolation is deliberate: it is what lets us measure the master and not the solvers. The inputs of each FMU are the other mass’s position and velocity.

The script implements Jacobi and Gauss-Seidel, optional linear extrapolation of inputs from the previous two outputs, and optional fixed-point iteration with rollback (state snapshot, then restore, then re-step with improved end-of-step input estimates, linearly interpolated across the step).

import numpy as np
from scipy.linalg import expm

m = (1.0, 2.0); k = (20.0, 30.0); kc, c = 40.0, 0.2
T = 10.0

def ref(t):                       # monolithic exact solution
    A = np.array([[0, 1, 0, 0],
                  [-(k[0]+kc)/m[0], -c/m[0], kc/m[0], c/m[0]],
                  [0, 0, 0, 1],
                  [kc/m[1], c/m[1], -(k[1]+kc)/m[1], -c/m[1]]])
    return expm(A*t) @ np.array([1.0, 0, 0, 0])

class Fmu:
    """x'' = (-(k+kc) x - c v + kc xo + c vo)/m ; input (xo,vo) linear over a step."""
    def __init__(s, i, x0):
        s.x = np.array([x0, 0.0]); mi, ki = m[i], k[i]
        s.A = np.array([[0, 1], [-(ki+kc)/mi, -c/mi]])
        s.B = np.array([[0, 0], [kc/mi, c/mi]])
    def do_step(s, h, u0, u1):    # exact for a linear input signal
        M = np.zeros((5, 5)); M[:2, :2] = s.A; M[:2, 2:4] = s.B
        M[2:4, 4] = (u1 - u0) / h
        s.x = (expm(M*h) @ np.concatenate([s.x, u0, [1.0]]))[:2]

def run(H, scheme, extrap=False, iters=0):
    f = [Fmu(0, 1.0), Fmu(1, 0.0)]; prev = None; err = 0.0
    for n in range(int(round(T/H))):
        o = [f[0].x.copy(), f[1].x.copy()]      # outputs at t_n
        snap = [o[0].copy(), o[1].copy()]       # getFMUState
        g = [o[1] + (o[1]-prev[1]) if (extrap and prev) else o[1],
             o[0] + (o[0]-prev[0]) if (extrap and prev) else o[0]]
        if scheme == "jacobi":
            f[0].do_step(H, o[1], g[0]); f[1].do_step(H, o[0], g[1])
        else:                                    # Gauss-Seidel: FMU 1 sees NEW output of FMU 0
            f[0].do_step(H, o[1], g[0]); f[1].do_step(H, o[0], f[0].x.copy())
        for _ in range(iters):                   # fixed point with rollback
            e = [f[0].x.copy(), f[1].x.copy()]
            f[0].x, f[1].x = snap[0].copy(), snap[1].copy()   # setFMUState
            f[0].do_step(H, o[1], e[1])
            f[1].do_step(H, o[0], e[0] if scheme == "jacobi" else f[0].x.copy())
        prev = o
        r = ref((n+1)*H)
        err = max(err, abs(f[0].x[0]-r[0]), abs(f[1].x[0]-r[2]))
    return err

I ran this with numpy 2.5.3 and scipy 1.18.1, measuring the maximum absolute position error of either mass at the communication points over 10 s, for five step sizes. Results (the numbers below are the actual script output, in metres; the motion amplitude is about 1):

H (s)	Jacobi, held inputs	Gauss-Seidel	Jacobi + linear extrap.	Gauss-Seidel + extrap.	Jacobi + 1 iteration	Jacobi + 3 iterations
0.1	5.13e+01	3.24e+00	3.53e+00	7.04e-01	2.38e-01	2.51e-01
0.05	3.82e+00	8.45e-01	3.59e-01	1.40e-01	6.92e-02	6.44e-02
0.025	8.55e-01	3.50e-01	8.52e-02	3.42e-02	1.70e-02	1.65e-02
0.0125	3.52e-01	1.71e-01	2.09e-02	8.50e-03	4.16e-03	4.14e-03
0.00625	1.73e-01	8.63e-02	5.20e-03	2.12e-03	1.03e-03	1.04e-03

Observed convergence order between the last two step sizes (log2 of the error ratio): 1.03 for Jacobi, 0.99 for Gauss-Seidel, 2.01 and 2.00 with linear extrapolation, 2.02 and 2.00 with one and three iterations. Three conclusions follow, with the usual caveat that this is one linear system with one parameter set, not a benchmark of any tool.

First, the plain Jacobi and Gauss-Seidel schemes are both first-order, so halving H halves the error; Gauss-Seidel buys a constant factor of about two, not a better order. Second, errors above 1 in a system whose amplitude is about 1 are not “inaccurate”, they are unstable: at H = 0.1 the Jacobi run amplifies. The fastest mode here has a period of about 0.71 s (the other is 1.54 s, from the eigenvalues of the monolithic system), so H = 0.1 is only about seven steps per period for a scheme that lags every coupling. Third, the biggest wins come from information quality: linear extrapolation raises the order to two and cuts the error by a factor of about 33 at H = 0.00625 (0.173 to 0.0052). One iteration with rollback gets an order-two result with a smaller constant (0.00103) at the cost of two do-step calls per FMU per macro step.

What the table does not show is as important. The extra iterations beyond one bring nothing here (3 iterations are no better than 1 and sometimes slightly worse at H = 0.1 and 0.00625) because linear interpolation of the coupling signal across the step caps the achievable order at two; iterating only converges to the answer of the interpolated scheme, not to the true solution. This is a general warning: iteration removes the lag, not the interpolation error. To go beyond second order you need higher-order polynomial input representation, such as the Hermite approach via output derivatives, along with the iteration.

If you run this against a real FMU pair with FMPy, the structure is the same. FMPy supports FMI 1.0, 2.0 and 3.0 per its README and can install from PyPI with pip install fmpy[complete]. Its high-level simulate_fmu drives a single FMU, so a multi-FMU master is written by instantiating each FMU through the lower-level classes and calling their do-step methods yourself, which is exactly the pattern in the script above. The earlier FMPy FMI 3.0 tutorial walks through that instantiation.

Packaging with SSP, Tooling, and Real-Time Operation

SSP: making the topology portable

A master needs a machine-readable description of the system, and SSP is the vendor-neutral answer. An SSP archive (extension .ssp) is a zip file containing a System Structure Description (SSD) with the components (FMUs or nested systems), their connectors, and the connections between them, plus parameter files (SSV for values, SSM for mappings) and the resources themselves. Separating structure from parameter values is the quiet power of the format: the same topology can be run with a nominal parameter set, a calibrated set from the physical asset, and a fault-injection set, with only the SSV swapped.

Per the SSP 2.0 announcement, the 2.0 release adds architecture exchange (so a topology can be shared before the models exist and refined incrementally), direct support for FMI 3.0 based features for more robust co-simulation, support for Modelica components and acausal connections, mapping of structured interface types to binary connectors, virtual ECU integration, embedding of extra content such as 3D models for twin use, and more efficient parameter updating through dynamic array sizes. The announcement did not name which tools implement which parts at release, so check each tool’s compliance notes instead of assuming.

Figure 4: SSP separates topology (SSD) from values (SSV). Different masters can load the same package, which makes cross-tool verification possible.

That cross-tool verification is an underrated benefit. If two independent masters load the same SSP and produce trajectories that agree to within the expected coupling error, you have evidence the result is a property of the system and not of one vendor’s scheduling quirk. If they disagree by more than the error model predicts, you have found either a defect or an unspecified behaviour such as a differing step order.

What SSP does not standardise is the algorithm: it can carry hints and, in 2.0, more structure for co-simulation, but the choice of Jacobi or Gauss-Seidel, of step size, and of extrapolation still lives in the master’s configuration. Record it next to the package, versioned together, or the result is not reproducible.

Tools: what each actually offers

The tooling landscape is changing and I limit myself here to what primary documentation states.

FMPy (CATIA-Systems on GitHub) supports FMI 1.0, 2.0 and 3.0, Co-Simulation and Model Exchange, runs on Windows, Linux and macOS, and offers a command line, a GUI and a web app. It is an excellent reference and prototyping master but is not a ready-made multi-FMU orchestrator, so a Jacobi or Gauss-Seidel loop like the one above is yours to write.
OMSimulator is described by OpenModelica as an FMI and SSP-based co-simulation environment; its user guide includes running SSP files, exporting SSV and SSM templates, and solver and fixed or variable step settings. I could not confirm from the pages I fetched which FMI versions it currently supports or the details of its coupling algorithm, so confirm against its documentation before relying on FMI 3.0 specific features. The repository lists v2.1.3 on the 2.1 maintenance branch at the time I looked.
INTO-CPS Maestro 2 is a JVM co-simulation orchestration engine with a domain-specific language (MaBL), a web UI and a CLI. Its documentation, at the time of reading, covers FMI 2.0 co-simulation, lists a fixed-step Jacobi algorithm with rollback and a variable-step algorithm without rollback, mentions iterative initialisation for strongly connected components, and lists slow-down to real time and digital twin capabilities among features in progress. The docs title carries an alpha label but link a 2.0.2 release; check the current status yourself.

For a digital twin that must ingest live data, the practical pattern is to embed one of these engines behind a service interface and drive it from the twin’s runtime, or to write a thin purpose-built master like the toy above when you need control over scheduling and data conditioning.

Real-time and hardware-in-the-loop

Everything above assumed simulated time runs as fast as the CPU allows. For real-time simulation and HIL (a physical controller or sensor connected to the simulated plant), the master must also meet a deadline: every macro step must complete within h seconds of wall time, every time. That constraint reshapes the algorithm choices.

Variable steps are usually off the table, since the communication period is dictated by the hardware’s I/O rate; rollback is also off the table because a physical device cannot be rolled back, so there is no “try again” once an output has been sent to a real actuator. That leaves fixed-step schemes with predictive input conditioning, which is why extrapolation is more important in HIL than in offline work, and why the accuracy of the first-order scheme is often unacceptable at the I/O rates available. Jacobi has the practical edge here: all FMUs can run on separate cores, and the worst-case step time is the slowest FMU plus communication, a bound you can measure. Gauss-Seidel serialises the chain, so the worst-case step time is the sum.

Latency adds an error that is easy to forget: a real device sees outputs delayed by the macro step plus I/O time, so a HIL coupling is at least as delayed as a Jacobi one. Treat the delay as part of the plant model and check the loop’s stability margin with the delay included. Where a step overruns, decide the policy up front: skip and hold the last output (safe for slow plants), run a shorter catch-up, or fail the test. Silently running late invalidates the result, and the master should count and report overruns.

Maestro’s feature list names slowing down to real time as a work in progress, which is a reminder that pacing against a wall clock is not free in a general orchestration engine: the sleep, jitter, and priority behaviour of the operating system matter, and a general-purpose OS with garbage collection (as on the JVM) is soft real time at best. Hard-real-time HIL typically runs on a real-time operating system or dedicated hardware, with FMUs compiled to native code.

Trade-offs, Gotchas, and What Goes Wrong

Stability is not accuracy. The classic failure is tuning the step size until the plots look right. A coupled mechanical system with a stiff coupling element can be unstable for a Jacobi step that an engineer would call small, and the failure is a slow growth of oscillation rather than a crash. Always compare against energy: if total mechanical energy of a conservative coupled system grows, the master is injecting energy.

Stiff coupling is the hard case. The toy system uses a moderate spring. Make the coupling much stiffer relative to the masses and the stable step size shrinks in proportion, because the coupling frequency sets the limit. At the extreme (a rigid connection) the coupling becomes an algebraic constraint, and explicit non-iterative co-simulation cannot handle it without a coupling formulation designed for it (for example, exchanging forces on one side and motion on the other, which is the usual force-displacement vs displacement-force choice).

Which variable crosses the interface matters. Passing positions both ways, as in the toy, is one causality choice. Exchanging force from A to B and position from B to A is another, and the stability regions differ. This is the reason a master cannot be separated from the modelling of the interface.

Event handling breaks naive masters. A discontinuity in an FMU inside a macro step can be missed or smeared. FMI 3.0’s early return and event mode exist for this, but the master has to use them. Without them, the effective event time is quantised to h.

Rollback assumptions hide state. Restoring an FMU’s state does not restore a random number generator the FMU shares with another, an external file handle, or a network socket. Deterministic replay needs FMUs that keep everything in the serialised state.

Parallel does not mean fast. Jacobi parallelism pays only if per-step compute outweighs synchronisation cost. For FMUs that step in microseconds, the thread hand-off dominates and a single-threaded Gauss-Seidel loop is faster in wall time.

Check the dependency declarations. A supplier FMU that declares no direct feedthrough but actually has some produces a hidden algebraic loop, which a Gauss-Seidel master will treat as a lag and which then shows up as an unexplained accuracy ceiling. Test with a step-size sweep like the table above: if the observed order is stuck at one when you expected two, something is misdeclared.

Unverified numbers are common. Published comparisons of master algorithms depend heavily on the benchmark system, and many rest on one case. Treat any single table, including mine, as an illustration of a mechanism, not a ranking.

Practical Recommendations

Start from the system’s topology, not from a favourite algorithm. Load the connection graph, read each FMU’s declared dependencies, and compute the strongly connected components. Everything outside a component can be solved by plain serial ordering with no iteration. Everything inside one needs a decision.

For offline analysis with few units and tight coupling, use Gauss-Seidel with linear extrapolation, and enable step-size control with rollback if every FMU supports state restore. For many loosely coupled units, or when the hardware offers cores to exploit, use Jacobi, but always with extrapolation, and size the step from the fastest coupled mode (a first sanity rule is at least twenty steps per period of the fastest coupled oscillation, a rule of thumb of mine rather than a standard). For real-time and HIL, fix the step, use Jacobi, extrapolate, and budget the worst-case step time with margin.

Measure your own coupling order. Run the system at three step sizes, compute the error against a fine-step run (or a monolithic reference if one exists), and fit the order. If it is one when you wanted two, you have a data-conditioning problem and not a solver problem.

Checklist:

[ ] Store the SSP, the master configuration (scheme, step policy, extrapolation order) and tool versions together.
[ ] Read declared input-output dependencies and find strongly connected components before running.
[ ] Confirm which FMUs support canGetAndSetFMUState before choosing a rollback-based design.
[ ] Run a three-point step-size sweep and report the observed order.
[ ] Monitor a conserved or bounded quantity (energy, mass) as an instability alarm.
[ ] For HIL, count deadline overruns and fail the run when a threshold is exceeded.
[ ] Cross-check at least once with a second master loading the same SSP.

Frequently Asked Questions

What is a co-simulation master algorithm?

A co-simulation master algorithm is the controller that advances several independently solved simulation units through time. At each communication point it passes outputs to inputs, calls each unit’s step function, and decides whether to accept, retry or resize the step. It determines how much coupling error accumulates, whether the run is stable, and how much parallelism is available. FMI standardises the unit interface but not the master.

Is Jacobi or Gauss-Seidel better for co-simulation?

Neither universally. Jacobi steps all units from the same snapshot, so it parallelises and has bounded step time, but every coupling lags by a step. Gauss-Seidel runs units in sequence so later units see fresh outputs, which typically gives lower error at equal step size, at the cost of serial execution and order sensitivity. In my two-mass test, both were first-order and Gauss-Seidel had about half the error of Jacobi.

Do I need rollback with getFMUState and setFMUState?

Only for variable-step control based on step doubling, for fixed-point iteration within a step, and for event location by bisection. Fixed-step masters with extrapolation need no rollback. Rollback needs every FMU to declare that it can get and set state, and many commercial FMUs do not. It also cannot undo external side effects, and it is impossible with physical hardware in the loop.

How do you handle algebraic loops in FMI co-simulation?

First, read the declared dependencies to confirm a loop really has direct feedthrough on every edge. Then choose: break it with a one-step delay and accept the error, iterate at a fixed time with fixed-point or Anderson-accelerated updates, or restructure the models so one exposes the implicit equation. Plain iteration converges only if the loop gain magnitude is below one, whereas Anderson acceleration does not require that.

What is SSP and how does it relate to FMI?

System Structure and Parameterization (SSP) is a Modelica Association standard that describes a composite system: components, connectors, connections and parameter sets, packaged in a portable archive. FMI describes a single model unit; SSP describes how units are wired together. SSP 2.0 was announced in late 2024 with architecture exchange and FMI 3.0 based features. It standardises topology, not the master algorithm.

Can FMI co-simulation run in real time for hardware-in-the-loop?

Yes, with constraints. The macro step must be fixed and every step must finish within its period, so variable-step control and rollback are unavailable. Jacobi stepping with input extrapolation across cores is usual, FMUs are compiled to native code, and the loop runs on a real-time operating system or dedicated hardware. Account for I/O delay as part of the plant, and count missed deadlines explicitly.

Co-Simulation Master Algorithms: Jacobi vs Gauss-Seidel, SSP and Algebraic Loops for Digital Twins (2026)