Variable Structure System Simulation via Predeﬁned Acausal Components

This article outlines a new approach of the experimental open-source modeling and simulation system Modia to simulate systems where the number of variables and equations can be changed after compilation and also during simulation, without having to re-generate and re-compile the code. Details are given for heat transfer in an insulated rod, where the discretisation of the rod is completely hidden from the symbolic engine. It is discussed how this approach could also be used in a future version of Modelica and/or FMI. Furthermore, this feature is also used in various variants to speed up collision handling in 3D mechanical systems. For example, by rigidly ﬁxing an object after it has been gripped, with or without calculating the elastic response, and thereby dynamically changing the number of degrees of freedom.


Introduction
Modia (Elmqvist et al. 2021) is an experimental, open source modeling and simulation system to develop new approaches to overcome the limitations of declarative, equation-based modeling languages such as Modelica (Modelica Association 2023).Modia is implemented with the powerful Julia programming language (Bezanson et al. 2017).It consists of a set of Julia packages, in particular of Modia.jl1 for equationbased modeling à la Modelica and of Modia3D.jl 2 for modeling of multibody systems.Neumayr and Otter (2023) extend Modia to process so called predefined acausal components 3 .These model components consist of a (usually small) set of Modia equations in which Julia functions are called that contain the core variables and equations of the components.These variables and equations can appear and disappear during simulation, without re-generation and re-compilation of code and without knowing in advance which model equations are utilized during such a simulation.
In contrast to this new approach, all previous proposals for systems with variable structure must either know in advance the entire models for all modes and switch between these models during simulation, (e.g., Mehlhase 2014;Mattsson, Otter, and Elmqvist 2015;Tinnerholm, Pop, and Sjölund 2022).Or the entire model is newly processed and code is re-generated and re-compiled (or interpreted) whenever the equation structure is changed4 , (e.g., Zimmer 2010;Tinnerholm, Pop, and Sjölund 2022).
In this article, the novel approach in Modia for modeling predefined acausal components is demonstrated with 1D heat transfer in a rod, where the number of discretization nodes can be changed before simulation start without re-compilation.It would also be possible to change the number of discretization nodes during simulation.
Modia3D is a more complex predefined acausal component.It was recently extended to cope with variable structure systems where the number of degrees of freedom can change during simulation, without re-compilation.The core part of this article discusses how this new feature is used to improve collision handling as an extension to the elastic response calculation introduced in (Neumayr and Otter 2019).Elmqvist et al. (2021, Section 2) describe the Modia Language, and Neumayr and Otter (2023, Appendix A) provide a short overview of it.A more detailed explanation is available in the Modia Tutorial5 .

Predefined Acausal Components
Neumayr and Otter (2023) introduce predefined acausal components which are based on a proposal of Elmqvist (2022): The equations of an acausal component are split into causal and acausal partitions.The intuition is that the causal partition is always evalu-ated in the same order, regardless of how the component is connected with other components.This partition is sorted, explicitly solved for the unknowns, and implemented with one or more functions.In contrast, sorting and solving of the acausal partition depends on the actual connection of the component.This partition is kept as a set of equations.Note that, the figures and some text fragments used below in this section are from Neumayr and Otter (2023) In Neumayr and Otter (2023), various variants of this basic approach are discussed.In particular, (a) the variables computed in the causal partitions can either be still visible in the equation part as proposed by Elmqvist (2022) or (b) a large part of these variables is hidden in the functions and do no longer appear in the equation part.Variant (a) has the advantage, that index reduction is still possible by differentiating the functions of the causal partitions.Variant (b) has the advantage that the causal partition, in particular the number of its variables and equations, can be changed after compilation and during simulation.Variant (b) has the drawback that index reduction is no longer possible for the causal partitions.Index reduction in the acausal partitions is still possible and is sufficient in many practical cases.However, it is not possible to use a predefined, acausal component as an inverse model if implemented with variant (b).In Modia and Modia3D variant (b) is used.
In Figure 1 the communication structure between the solver, the sorted and solved equations and the functions6 of the causal partitions are shown: The variables of the solver (state vector x and the vector of event indicators z) are split into an invariant and a variant part: x = (x inv , x var ) and z = (z inv , z var ).The dimensions of the invariant parts are fixed before the simulation starts.The dimensions of the variant parts, which are contained in the functions of the causal partitions, can change at events during simulation.Since x var , z var are communicated directly between the functions and the solver, the symbolic processing of the equation part of a model is not affected by these variables.Therefore, these variables can in principle be changed at event times -variables can be added or can disappear.
This basic approach is demonstrated using the predefined acausal component of Figure 2, which models heat transfer in a rod with an insulated surface.On the left and right sides of the rod, thermal connectors a, b are present (called port_a, port_b in Listing 1) with potential variables a T , b T (temperatures) and flow variables a Q flow , b Q flow (heat flow rates).The partial differential equation, which mathematically describes the heat transfer in one dimension is discretized in space by volumes V i = ∆x • A of equal sorted and solved equations functions of predefined acausal components Communication between the solver, the sorted and solved equations, and the functions of the predefined acausal components.The state vector x and the event indicators z are split into an invariant and a variant part: x = (x inv , x var ), z = (z inv , z var ).The variant parts consist of the states defined and used in the causal partitions of all predefined acausal components.The dimensions of the invariant parts are fixed before simulation starts.The dimensions of the variant parts can change at events during simulation.lengths ∆x and identical areas A. In the center of volume i, a temperature T i is defined, leading to a temperature vector T = [T 1 , T 2 , . . ., T n ].
Listing 1. Simple usage of insulated rod InsulatedRod2 with one-dimensional heat-transfer.On the left side it is connected with a fixed temperature source FixedHeatFlow with T = 220 °C = 493.15K, and on the right side with a fixed heat flow source FixedHeatFlow with Q flow = 0. # plot temperatures plot ( heatedRod, ( " fixedT.port.T " , " rod.T " ) ) In Listing 1, a Modia model is shown with a predefined acausal component InsulatedRod2 of the rod7 .Its left thermal connector port_a has a fixed temperature source FixedTemperature.Its right thermal connector port_b has a fixed heat-flow source Fixed-HeatFlow with the default zero heat-flow rate.This means that, the rod is completely insulated on the right side and has a fixed temperature on the left side.Note that, A|B merges model or parameters B with model A. Command @instantiateModel(Heated-Rod) symbolically processes this model and generates Julia code that is translated to executable code.The simulate!statement changes the discretization, and thus the dimension of the temperature vector T , from 5 to 8 volumes before simulation starts without a The implementation of the InsulatedRod2 model is shown in Listing 2. To start with, file Insulated-Rod2.jl is included containing the definition of a Julia struct holding the data and the local variables of the component, as well as some Julia functions.More details are given below.A standard Modia definition of the model is then given defining the parameters and connectors of the component.Contrary, to a standard Modia component, no equations are present.For simplicity, no units are used in this model and its associated functions.However, the actual implementation of the component in Modia supports units.The component is a predefined acausal component model, since the following special parameters are provided at the beginning of the model, defining functions to be used for symbolic processing and during simulation: • _buildFunction: Before symbolic processing begins, the hierarchical dictionary of the root model to be compiled is traversed and function buildInsulatedRod2, defined from _build-Function, is executed for each sub-model it contains.This function (a) defines additional model variables and equations that are merged with the corresponding model and (b) returns an instance of the Julia structure, which acts as the internal memory of the component.
• In Listing 3, the implementation of function build-InsulatedRod2 is shown.In this function, the model instance (of the actual InsulatedRod2 component) is merged with additional model definitions consisting of two variables and several equations.It returns the merged model.Additionally, the internal memory of the component is instantiated with Insulated-RodStruct() and is also returned by buildInsulated-Rod2!.This internal memory is later identified by the unique identifier ID, which is specified in the function call.Function call openInsulatedRod! in the equation section copies the rod temperatures T from the state vector of the simulation engine into the InsulatedRod-Struct memory and returns a reference to it as ins-Rod.The argument list of this function call includes the unique identification ID of the predefined acausal component.It is provided when buildInsulatedRod2! is called.$ID is inside an Abstract Syntax Tree, due to :[...] and $ inserts the actual (literal) at this place.In Julia terminology this is called "interpolation".The insRod reference is then used in subsequent function calls, for example, to retrieve the value of the first temperature node with getT1(insRod).This value is used in an equation to calculate the heat flow from port_a to the first internal node.Finally, com-puteInsulatedRodDerivatives! computes the derivatives of the temperatures and copies them into the state derivative vector of the simulation engine.As can be seen, the equation section is independent from the number of discretization elements nT.Therefore, the number of these discretization elements can be changed without re-generation and re-compilation.The implementation of the _initSegmentFunction is shown in Listing 4. This function is called before a new simulation segment is initialized.The first statement inquires the reference insRod of the internal memory of the component.Before the first simulation segment, the (evaluated) parameters are stored in insRod.Furthermore, some dependent parameters are computed and also stored in this memory.Finally, new_x_segmented_variable is called to define the name of a new state, its derivative and its unit together with the initial value of insRod.T which is the current value of vector T. It is initialized with parameter T_init in initFirstSegmentInsulatedRod2!.Note that, even if the InsulatedRod2 component always uses the same definition, the states must be newly defined for each new simulation segment.Note that, a similar approach could be implemented in Modelica with reasonable effort: The simplest implementation would be to use External Objects and add additional utility functions (Modelica Association 2023, Section 12.9.6-12.9.7).These are equivalent to the utility functions of Modia, e.g., to add variables at events.These utility functions would directly communicate with the underlying simulation engine.If they are available, the Modia example of the insulated rod could be implemented as outlined in Listing 6.The main benefit would be that the number of temperature nodes can be changed after translation and that the Modelica model consists essentially of three scalar equations.These equations are independent from the number of temperature nodes.The drawback is that functions openInsulatedRod, close-InsulatedRod, getGe2, getT1, getTend, computeInsu-latedRodDerivatives need to be implemented in C. Note that these would be simple C-functions.For example, computeInsulatedRodDerivatives could be implemented as shown in Listing 7.
) ) ; } c o p y _ d e r _ x _ s e g m e n t e d _ v a l u e _ t o _ s t a t e ( m, insRod-> T_startIndex, der_T ) ; return 0; } The eFMI standard (Functional Mockup Interface for embedded systems, (Lenord et al. 2021)) defines an intermediate language GALEC to transform acausal models to production code.GALEC is basically a very small subset of the Modelica language with some extensions as needed for embedded systems.The extension also includes a simple form of member functions.If such member functions were supported in Modelica, the implementation of predefined acausal components could be done completely in the Modelica language and without External Objects or C-code.
Modelica models can be exported in FMI format (Modelica Association 2022).This includes Modelica models with External Objects.The FMI standard communicates the values of variables explicitly with setter and getter function calls.In principle, it would be possible to add another variable type to the FMI standard, where internal variables (including states) are communicated directly to the solver and no longer via the setter/getter function calls.If such an enhancement were available, the causal partition part of the Modelica model of Listing 6 could be transformed to an FMU, where the node temperatures would no longer be communicated via the FMI setter/getter functions and would no longer be visible in the environment in which the FMU is used.This would have the great advantage that the number of variables and equations can be changed during the simulation.

Segmented Simulation and Collision Handling
Modia3D is a multibody tool for 3D mechanical systems implemented as a predefined acausal component of Modia according to section 2. Modia3D is targeted for solvers with adaptive step size control to compute results close to real physics including collision handling using the Minkowski Portal Refinement (MPR) algorithm (Snethen 2008;Neumayr and Otter 2017) and collision response for elastic contacts (Hertz 1896;Flores et al. 2011;Neumayr and Otter 2019).Modia3D has a very flexible and modular design pattern.It is extended (since v0.12.0) to cope with variable structure systems where the number of degrees of freedom (DoF) can change during simulation, without re-compilation.Modia3D offers two kinds of joints: The first kind of joints contains Modia equation sections with invariant variables, including invariant states.These invariant elements are visible for Modia and cannot be removed or added during simulation.The interface to the Modia3D functionality is designed to define differential equations only on the Modia side in Modia equation sections, so that state constraints can be defined and index reduction can be performed on invariant states.The joints of the second kind define variant variables, including variant states, which are visible only in the Modia3D predefined acausal component.These joints can be added or removed during simulation.For example, an Object3D has an optional keyword fixedToParent with a default value of true.In this case, the Object3D is rigidly connected to its parent Object3D.This means it has zero degrees of freedom.If the value is set to false, the Object3D is allowed to move freely with respect to its parent, meaning it has 6 degrees of freedom and 12 variant states.At events, keyword fixedToParent can be changed from false to true and vice versa.Neumayr and Otter (2023, Table 2) define Modia3D actions which modify the second (variant) kind of joints and trigger structural changes during simulation, e.g., actionAttach, actionReleaseAndAttach, actionRelease, actionDelete.The new states (joints) added during simulation with e.g., actionRelease are initialized based on the last known position, velocity, acceleration and rotation.All remaining states are re-initialized with their last known values.Based on that, the internal 3D structure is rebuilt and executed until another action for a structural change is triggered.This restructuring is performed with dynamic data structures and is extremely fast (< 1 ms).If the sphere is rigidly attached to the plate or gripper, there are no states, and nothing is displayed.The sphere in Scenario 4 (S4) has no states.Therefore, nothing is displayed.The states for scenarios 2 and 3 (S2, S3) are available and are displayed if the sphere is freely moving.For the absolute position of the sphere center see Figure 6.In this section, several combinations of segmented simulation and collision handling are discussed using a KUKA YouBot robot gripping and transporting a sphere, see Figure 4.This robot has a 5 DoF arm and was manufactured in the years 2010-2016.Elmqvist et al. (2021) model the drive trains and controllers of the robot in Modia, and the 3D mechanics with Modia3D.
Four variants of the following transportation scenario are simulated.In all these scenarios, the robot follows the same trajectory.Initially, the cargo, e.g., a sphere, rests on a plate.It is gripped by the robot's gripper and transported upwards until it is placed down again, where it rests on the plate until it is 8 Intel(R) Core(TM) i7-9850H CPU @ 2.6 GHz, RAM 32 GB gripped again.The states of the freely moving sphere (see Figure 5), if available, and the absolute position of the sphere center (see Figure 6) are displayed.The simulation times of all four scenarios are compared in Table 1.
Scenario 1 (S1)9 : The transportation scenario is modeled with collision handling, compare (Neumayr and Otter 2023, Scenario 2(b)).This means, the sphere collides with the plate, as well as with the fingers of the gripper.
Scenario 2 (S2)10 : The transportation scenario is modeled with segmented simulation and collision handling, see Listing 8. DoFs are added or removed during simulation: At the beginning, the sphere is rigidly attached to the plate.Shortly before the gripper reaches the sphere, the sphere is released (+6 DoF) and collides with the plate.Shortly afterwards it collides with the gripper.After approximately one second, the sphere is rigidly attached to the gripper (-6 DoF).Until the gripper is again close to the plate to release the sphere (+6 DoF), which collides with the plate.Collision handling remains on even if the sphere is rigidly connected to the gripper or plate, as collisions with other bodies can still occur.Scenario 3 (S3)11 : The transportation scenario is modeled with segmented simulation and collision handling.Scenario 3 is very similar to Scenario 2, except that collision handling is disabled when the sphere is rigidly connected to the gripper or plate, since in this scenario it is already known that no further collisions will occur.This is deactivated with enableContact-Detection = false in Listing 8. Basically, this means that the distance calculations between each collision pair is switched off in these phases.
Scenario 4 (S4)12 : The transportation scenario is modeled with segmented simulation only, compare (Neumayr and Otter 2023, Scenario 2(a)) and Listing 9. Collision handling is switched off for this scenario.This means, the sphere is rigidly attached to the plate, when resting, and rigidly attached to the gripper during transportation.Each time the sphere is rigidly connected to the plate or gripper, the segment is re-initialized.Since the relative velocity and angular velocity between the sphere and the gripper is zero, when the sphere is attached to the gripper or attached to the plate, the physics is correctly modeled under the idealized assumption that gripping time is infinitely small.Basically, this means that gripping effects are neglected.The simulation time of Scenario 4 is about 19 times less than that of Scenario 1.This is because Scenario 4 (segmented simulation only) is basically a non-stiff system where the solver can use large step sizes.In addition, the time for reconfiguration of the multibody system, for gripping and releasing, is negligible.Fine-tuning of collision handling during transportation of the gripped freight is no longer required.Furthermore, any type of cargo can be transported, regardless of its shape.The disadvantage is that the details of the gripping are not modeled, but this can be important.
Scenario 1 (collision handling only) is a stiff system because the gripper holds the sphere by elastic contact and friction forces, which change during transport.Therefore, the solver must use much smaller step sizes.One limitation of collision handling with the MPR algorithm is that it only supports point contacts.If the cargo would be a box, see (Neumayr and Otter 2023, Scenario 3(b)), it would not be possible to calculate a unique point contact that is continuous over time, for example, because one box and one gripper face or one box and one plate face are parallel to each other during contact.All these considerations lead to a compromise in modeling the gripping and releasing of the cargo with collision handling, and otherwise rigidly attaching the sphere to the plate or gripper, resulting in Scenario 2 and Scenario 3.
There is not such a big difference in simulation time for Scenarios 1,2,3, see Table 1.In all three cases, the calculation of the elastic contact response is the limiting factor.This effect is modeled in all these cases.In more realistic scenarios, the approach of Scenario 2 or 3 may pay off, if the number of collision phases is small relative to the remaining actions.

Conclusion
The novel approach of variable structure systems with Modia/Modia3D seems to be very promising.De- pending of the predefined acausal component and the application one can design (extend) specific actions to trigger new segments and re-initialize the model.In this paper, existing Modia3D actions are extended by enabling or disabling collision handling during simulation, which speeds up the simulation and allows to model form locked fixing of cargos.Furthermore, an example was sketched of how the Modelica language and the FMI standard could be enhanced to allow the number of variables and equations to be changed during simulation.

Figure 2 .
Figure 2. Space discretized partial differential equation of one-dimensional heat transfer in a rod with an insulated surface.It is defined with parameters L (length of rod), n (number of volumes), A (area), (density), c (specific heat capacity), λ (thermal conductivity), T 0 (initial value in each volume), states T i (temperatures in the center of each volume), thermal connectors a, b with potential variables a T , b T (temperatures), and flow variables a Q flow , b Q flow (heat flow rates).

Figure 3 .
Figure 3. Plot of temperatures of heated rod model.

Listing 4 .
_initSegmentFunction definition.# Called once before new sim.segment function i n i t S e g me n t I n s u l a t e d R o d2 !( m, path, ID, parameters ) insRod = g e t _ i n st a n t ia t e d S ub m o d e l ( m,ID ) if i s F i r s t I n i t i a l O f A l l S e g m e n t s ( m ) i n i t F i r s t S e g m e n t I n s u l a t e d R o d 2 !( insRod ; parameters... ) end # Define new states and state derivat.ins Ro d.T_ sta r tIn de x = n e w _ x _ s eg m e n t e d _ v a r i a b le !( m, path * " .T " , path * " .der ( T ) " , insRod.T, " K " ) return nothing end

Listing 7 .
C-function to compute state derivatives int c o m p u t e I n s u l a t e d R o d D e r i v a t i v e s ( struct M * m, struct InsRod * insRod, double Ta, double Tb ) { double * T = insRod-> T ; double * der_T = insRod-> der_T ; double

Figure 4 .
Figure 4. YouBot gripping or releasing a sphere on a plate.

Figure 5 .
Figure5.States of the sphere.They are equivalent to the translation of the sphere center in x, y, z direction with respect to its parent.If the sphere is freely moving, world is its parent.States can only be displayed if they are present.If the sphere is rigidly attached to the plate or gripper, there are no states, and nothing is displayed.The sphere in Scenario 4 (S4) has no states.Therefore, nothing is displayed.The states for scenarios 2 and 3 (S2, S3) are available and are displayed if the sphere is freely moving.For the absolute position of the sphere center see Figure6.

Figure 6 .
Figure 6.Absolute position of sphere center for the 4 scenarios.
FullRestart event before the root model is re-initialized at a new simulation segment.In both cases, all local variables of the predefined acausal component model (including states and zero-crossing functions) must be redefined, as well as initial values for newly defined states.
_initSegmentFunction: This function is called by the simulation engine before the root model is initialized and at each y _ d e r _ x _ s e g m e n t e d _ v a l u e _ t o _ s t a t e ( m, i n s R o d .T _s tar tIn dex , insRod.der_T ) return true endIn Listing 5, the most important remaining functions are shown.The state derivatives are computed and copied to the state derivative vector of the simulation engine with computeInsulatedRodDerivatives!.

Table 1 .
Mean x and standard deviation s of the simulation time of all four scenarios (S1-S4) each for n = 12 runs on a standard notebook 8 .
Listing 8. Robot program of Scenario 2. Collision handling is enabled by default.It can be turned off or on again in all action commands with enableContactDetection.Scenario 3 is defined, by setting this flag to false, as indicated in the comment lines.
Listing 9. Robot program of Scenario 4. Ac ti on R el ea s e A n d A t t a c h ( actions, " sphereLock " , " r ob o t. b as e .pl at e Loc k " )