S-Function Library for Bond Graph Modeling

S-functions are short for system-functions. They are used for extending the capabilities of Simulink. S-functions allows us to add our own algorithms to Simulink models. The process of creating S-function blocks is quite simple. Simulink provides a S-function API which can be used to write a S-function routine observing a set of laid down rules. The compiled routine is enclosed inside a Simulink block, which can be subsequently customised by masking. A library of customised S-function blocks are created for an application specific task. This library can be subsequently distributed to work in MATLAB environment.

If the system model has been already modeled as a set of mathematical equations, it becomes easy to convert each equation into a S-function block, and and develop the system model in Simulink.

Vectors in S-function
In a Simulink block, a vector of input, u, are processed by a vector of states, x to output a vector of output, y ( Fig. 1) (Mathworks, 2011). The state vector may consist of continuous states, discrete states, or a combination of both. u (input) x (states) y (output)

Fig. 1. Generalised Simulink block
In S-functions written in the MATLAB programming language, the MATLAB S-function, Simulink partitions the state vector into two spaces. The first part of the state vector is occupied by the continuous state, and the discrete states occupy the second part. But in the other programming language written S-function, MEX-file S-functions, there are two separate state vectors for the continuous and discrete states.

Steps in simulation
Routines have to be written in S-function to carry out the simulation steps required by the Simulink engine. To cross-reference the routine required for the simulation step, two different approaches are used. For an MATLAB S-function, Simulink passes a flag parameter to the S-function. The flag indicates the current simulation stage. Routines in M-code calls the appropriate functions for each flag value. For a C MEX S-function, Simulink calls the S-function routines directly. This is done by following a naming convention for the routines.

Steps in S-function block simulation
Simulink engine first calls the S-function Routine to perform initialisation of all S-functions block in the model. Later, the engine makes repeated calls during simulation loop to each S-function block in the model, directing it to perform tasks such as computing its outputs, updating its discrete states, or computing its derivatives. Finally, the engine invokes a call to each S-function Routine for a termination task (Fig. 2). The tasks required at various stages, include: • Initialization: Simulink initializes the S-function as a first step. The tasks are: -Initialising the SimStruct, a structure that contains information about the S-function.
-Setting the number and size of input and output ports.
-Setting the block sample times.
-And allocating storage areas and the sizes array. • Calculation of next sample hit -for a variable step integration routine, this stage calculates the time of the next variable hit, that is, it calculates the next stepsize in the variable step. • Calculation of outputs in the major time step. After this call is complete, all the output ports of the blocks are valid for the current time step.
• Update discrete states in the major time step. In this call, all blocks should perform once-per-time-step activities such as updating discrete states for next time around the simulation loop. • Integration: This applies to models with continuous states and/or nonsampled zero crossings. If S-function has continuous states, Simulink calculates the derivative of the continuous state at minor time steps. Simulink computes the states for S-function. If S-function (C MEX only) has nonsampled zero crossings, then Simulink will call the output and zero crossings portion of S-function at minor time steps, so that it can locate the zero crossings.

Flags in MATLAB S-function
An MATLAB file that defines an S-Function block must provide information about the model. Simulink needs this information during simulation. The MATLAB S-function has to be of the following form: [sys, x0, str, ts]= f (t, x, u, flag, p1, p2, ...) where f is the name of the S-function. Simulink passes t, the current time, x, state vector, u, input vector, integer flag in argument to S-function. In an MATLAB S-function flags are used for indicating the current simulation stage. S-function code calls the appropriate functions for each flag value. Table 1 lists the simulation stages, the corresponding S-function routines, and the associated flag value for MATLAB S-functions. An MATLAB S-function returns an output vector having the following elements: • sys -the values returned depend on the flag value (for flag = 3, sys contains the S-function outputs). • x0 -the initial state values at flag = 0 (otherwise ignored). • str -reserved for future use. • ts -two-column matrix containing sample time and offset.

C MEX S-function callback methods
As with MATLAB S-functions, Simulink interacts with a C MEX-file S-function by invoking callback methods that the S-function implements. C MEX-file S-functions have the same structure and perform the same functions as MATLAB S-functions. In addition, C MEX S-functions provide more functionality than MATLAB S-functions. C MEX-files can access and modify the data structure that Simulink uses internally to store information about the S-function. This gives it an ability to handle matrix signals and multiple data types.
C MEX-file that defines an S-Function block provides information about the model to Simulink during the simulation. Unlike MATLAB S-functions, no explicit flag parameter is associated with C MEX S-function routine. But the routines have to follow the naming convention. Simulink then automatically calls each S-function routine at the appropriate time during its interaction with the S-function. It defines specific tasks which include defining initial conditions and block characteristics, and computing derivatives, discrete states, and outputs.
Simulink defines in a general way the task of each callback and their sequence (Fig.3) (Mathworks, 2011). The green box in Fig.3 are the compulsorily present routines and the blue box are the optional routines. The S-function is free to perform any task according to the functionality it implements. For example, Simulink specifies that the S-function's mdlOutput method must compute that block's outputs at the current simulation time. It does not specify what those outputs must be. The callback-based API allows to create S-functions, and hence custom blocks, of any desired functionality. The contents of the routines can be as complex and any logic can reside in the S-function routines as long as the routines conform to their required formats.

MEX versus MATLAB S-functions
Both the approaches of MATLAB and MEX S-functions development have some inherent advantages due to their origin in the programming language they are coded in. The advantage of MATLAB S-functions is speed of its development. Developing MATLAB S-functions avoids the time consuming compile-link-execute cycle required when developing in a compiled language like C. MATLAB S-functions also have easier access to MATLAB toolbox functions and can utilize the MATLAB Editor/Debugger. MEX S-functions are more appropriate for integrating legacy code into a Simulink model. For more complicated systems, MEX S-functions may simulate faster than MATLAB S-functions because the MATLAB S-function has to call the MATLAB interpreter for every callback routine.

S-function examples
The ease of writing an S-function is demonstrated here by taking a simple example. Both MATLAB file and C Mex approaches will be demonstrated. A block 'timestwo' is implemented in S-function. This simple block has no states. Functionally, the block takes an input scalar signal, doubles it and outputs to a connected device (Fig.4).

MATLAB S-function implementation of timestwo
The S-function will be defined with four input arguments from Simulink. They are, the current time t, state x, input u, and flag flag. The routines that are called by timestwo.m are mdlInitializeSizes and mdlOutputs. The routine mdlInitializeSizes passes on the block information to the size structure. The information on number of outputs, inputs, continuous and discrete states, sample times and whether direct feed through is present, is passed on to the structure variable. The output variables, x0, str and ts are also set to the desired values. The second routine mdlOutputs just doubles the input scalar. The timestwo MATLAB S-function can now be used in the Simulink model, by first dragging an S-Function block from the User-Defined Functions block library into the model. Then entering the name timestwo in the S-function name field of the S-Function block's Block Parameters dialog box (Fig. 5).

C MEX S-function implementation of timestwo
C MEX S-function will contain the callback methods mdlInitializeSizes, mdlInitializeSampleTimes, mdlOutputs and mdlTerminate. The simulation steps will have setting of intial condition by the first two blocks, the simulation loop in the third block and final termination by the last block (Fig.6 ). The code is reproduced below: #ifdef MATLAB_MEX_FILE/ * Is this file being compiled as a MEX-file? * / #include "simulink.c" / * MEX-file interface mechanism * / #else #include "cg_sfun.h" / * Code generation registration function * / #endif

Explanation of C MEX timestwo code
We start by two define statements, and give the name to our S-function and tell the Simulink engine that the code is written in Level 2 format. Later we include the header file simstruc.h, which is a header file that defines a data structure, called the SimStruct, that the Simulink engine uses to maintain information about the S-function. A more complex S-function will include more header files.
The callback routine mdlInitializeSizes tells the Simulink engine that the function has no parameters, has one input and output port. It also declares that only one sample time will be specified later in mdlInitializeSampleTimes. Simulink is also informed that the code is exception free code. This declaration speeds up the execution time. The next callback routine mdlInitializeSampleTimes declares the sample time to be inherited i.e the block executes whenever the driving block executes.
The callback routine mdlOutputs calculates outputs at each time step. The input and output port signals is accessed through a vector of pointers. The width of the output port which is defined to be dynamically set is then read and the program loops for the width while calculating the output signal to be two times the input signal. The final callback mandatory routine mdlTerminate performs the end of the simulation task, which in the present case is a NIL set. following this routine the mandatory trailer code for compiler is present. Its absence will lead to compile errors.
The C MEX S-function has to be now compiled. Simulink gives a choice of C compiler to be used. We can use either the built in MEX compiler or any other C compiler already loaded in the system. The following command at the command line mex -setup allows us to locate all the compilers available in the system and the option to use one for compiling and linking in the MATLAB environment. Later the command mex timestwo.c compiles and links the timestwo.c file to create a dynamically loadable executable for the Simulink software to use. The resulting executable is referred to as a MEX S-function. The MEX file extension varies from platform to platform. For example, on a 32-bit Microsoft Windows system, the MEX file extension is .mexw32. The compiled run time file is then put into the S-function block similer to MATLAB S-function file (Fig. 5).

Analogous behaviour of physical systems
The bond graph approach to physical system modeling was conceptualized by Hank Paynter on April 24, 1959 (Paynter & Briggs, 1961), inspired by the earlier work of Gabriel Kron (Kron, 1962). Bond graph language is a port based graphical approach for modeling energy exchange between subsystems. This technique was further developed by Karnopp and Rosenberg (Karnopp et al., 1990;Karnopp & Rosenberg, 1968;Rosenberg & Karnopp, 1983). Several books, special issues and articles on bond graph technique have popularised it for growing usage (Borutzky, 2009;Borutzky et al., 2004;Breedveld, 1984;1991;2004;Breedveld et al., 1991;Cellier et al., 1995;Dauphin-Tanguy, 2000;Gawthrop, 1995;Gawthrop & Smith, 1996;Mukherjee & Karmakar, 2000;Thoma, 1990;Thoma & Perelson, 1976).  Behaviour of a physical system is constrained, either implicitly or explicitly by laws of physics viz. mass and energy conservation, laws of momentum and positive entropy production. Furthermore, various physical domains are each characterized by a particular quantity that is conserved. Each of these domains have analogous ideal behaviour with respect to energy (Table 2). This analogy led to the concept of energy port, the building block of bond graph modeling language. Here, the interaction between physical systems is through energy port and is always bidirectional. There will be an input signal and a consequent output signal ('back effect') and their product will signify the 'power that is transacted'. From a computational point of view, the effort could be computed by 'Port 1', while the flow is computed in 'Port 2'. It could be the other way around as well. Apriori the computational direction of signal is not known, except the fact that they are in opposite direction (Fig.7).

Bond graph elements
Bond graphs are labelled, directed graphs. The vertices of a bond graph denote subsystems, system components or elements, while the edges, called power bonds or bonds for short, represent energy flows between them. The nodes of a bond graph have power ports where energy can enter or exit. As energy can flow back and forth between power ports of different nodes, a half arrow is added to each bond indicating a reference direction of the energy flow. The amount of power, P(t), at each given time, t, is given by the product of the two conjugate variables, which are called effort, e, and flow, f, respectively.
There can be five groups of physical behaviour by elements handling energy: 1. Storing of energy. 2. Supply on demand. These five behaviour are represented by nine basic elements.
• Two types of storage elements, effort or Capacitive-storage and flow or Inductive -storage.
• Two types of sources, Sourceeffort and Sourceflow.
• Two types of Reversible transformators, non-mixing, reciprocal Transformer or TF-type transducer and mixing anti-reciprocal Gyrator or GY-type transducer. • Irreversible transducer is an energy Dissipater or can be also classified as entropy producing R-type transducer. • Distributor junction are also in dual form, 0-junction and 1-junction. The 0-junction represents a generalised domain independent Kirchoff current law and similarly a 1-junction represents a generalised domain independent Kirchoff voltage law.
The elements can also be segregated based on their port structure: • Five one-port elements: C, L(orI), R, S e , S f . • Two two-port elements: TF, GY.

Constitutive relations
A constitutive relation with a constant parameter characterises each element. For sources, the imposed variable is independent of the conjugate variable, and for the rest of elements, the relationship is algebraic between its conjugate variables. The storage elements are classified as memory elements. The preferred constitutive equation is integration with respect to time.
If differentiation with respect to time is used, the information on initial condition or history is lost. If sensors are included in the bond graph modeling and used for determining the factor of the algebraic relationship between the conjugate variable, we rename the bond as modulated bond with prefix M for the modulated element name.
Inductor, I    Table 4. Two-port elements

Computational causality
Each bond connects two power ports of different primitive elements and carries two power variables as can be seen in Fig.7. One of the two power variables may be determined by one of the two sub-models, while the other is determined by the other model. A short stroke, called causal stroke, perpendicular to the bond is placed at one of its ends of the bond. This indicates the computational direction of the effort variable.
Consequently the other open end is the decider of the flow variable. The nine basic elements with their constitutive relationships that are dependent on their causal stroke are shown in Table 3

S-function implementation for C bond
We will now use the C MEX S-function code to develop a continuous state Simulink block. An One-port element which stores energy is chosen for illustration. Effort causal Capacitor C element in Table 3

Junction algorithm
In a bond graph model a set of elements is connected at the junction. One of the elements in the set takes up the role of decider bond. Remaining bonds in the set per-se take up the role of non-decider bonds. The junction algorithm is illustrated by taking O-junction as an example. The 0-junction block is a common effort junction (Fig.8). The effort is decided by a decider bond attached to it and having the causal bar towards the junction. Similarly for a 1-junction, decider bond will have its causal bar away from the junction, thus complementing the 0-junction behaviour. In the figure, J in 's are the input from the non-decider bonds into the junction, and J sum is the output of the junction.
The governing law of the 0-junction, Effort Law (or KCL), states that the flow of the decider bond is the sum of the flows of all the non-decider bonds. In Fig.8, the decider bond of the junction has J sum , the sum of all J in , as its causal (flow) variable. Value of the conjugate variable (effort) of the decider bond, J out is decided by the bond's constitutive equation. The second part he governing law of the 0-junction, states that the all the bonds connected to the junction share a common effort. Thus the effort of the decider bond becomes the causal variable for all the non-decider bonds. For each non-decider bond, its non-causal variable  Table 6 tabulates the algorithm for one-port elements. Table 6. One-port element variables It can be seen that for an element connected to a junction, the causal variable can take either of the two values, J sum or J out , depending on whether the connection bond is a decider or a non-decider, respectively. Similarly the non-causal variable can have either the value of J out or J in , again depending on whether the connection bond is a decider or a non-decider.
A two-port is connected to two junctions at either end. These junction could be either 0-junction or 1-junction. But as a two-port element allows only two combinations of causality out of available four, there will be only 2 × 2 alternatives for a variable. Looking at the two-port junction in Fig.9, if the element is decider bond for junction J 1 , and the junction is 0-junction, then the flow at input port will be J 1sum . And if junction is a 1-junction the effort will be J 1sum . Similarly for the output junction connected to J 2 . Note that the direction of two-port half arrows is reversed in the two junctions. Thus if two-port element with flow causal at input port is decider bond for the first junction, the junction has to necessarily be an 0-junction, but flow casual at output port as decider bond indicates a 1-junction. The modulus of the two-port element decides the value of the variable on the conjugate port, while the junction decides the conjugate variable value of the port. In a similar manner for all combination of causality and decider bond, the variables can be listed out. The algorithm is summarised in Table 7 (Umarikar et al., 2006). Decider Non-decider Decider Non-decider Table 7. two-port element variables

Linking elements
A bond graph model on paper does not explicitly use a connector as in block diagram model, to link one element to another element. To retain the look and feel of a paper model when transferred to computer terminal, the connection between the elements has to be invisible. The masking properties of the Simulink block is utilised for this purpose in the tool box. A junction is considered as a node, to which the elements are connected. The bond graph element, as in paper model, is placed next to its associated junction. The link to the junction is made by entering the junction label in the mask parameter box. Shared Memory' algorithm is then used to implicitly connect the element.
In bond graph modeling two or more elements will be linked to a junction. Their data have to be shared. Shared data structure is used in the toolbox. The memory locations are earmarked for a junction by assigning it a unique label by character aggregation' during its first run. An S-function's initialisation callback method is used for memory allocation as this callback is used only once during the simulation run. The associated elements using the notation listed out in Table 6 and Table 7, share their respective memory address, thereby their data.
The elements in the tool box are masked and have screen interfaces. For a one-port element the following details need to be entered. For a two-port element the extra information of the second port is also entered. Similarly a junction screen interface will have all entries of the elements that are linked to it, along with their energy flow direction signs.

Propagation of data
The propagation of data from one element to the next is by reading and writing into a common block address (Fig.10). The input element, designated as Provider, produces the data through its constitutive equation and writes it on to the labelled memory address. The label, as discussed earlier is unique to a junction. The next element in hierarchy, designated as Consumer, in turn reads the data and by using its constitutive equations, Produces the next set of data which is written on to the next labelled memory address. There can be many Consumer of the data but there is only one unique Provider. An analogy to bond graph junction concept of one decider and many non-decider can be clearly seen here. It is also seen that a Consumer element of previous step becomes the Provider element in the next step. As the model is hierarchical, all the elements are in turn Provider and Consumer to their respective memory address block in one integration step. The memory location is released and freed when the data is not needed.

S-function implementation for shared memory link
The port of the bond graph element as discussed above is the pointer to the shared memory address. When we specify a input/output port for a bond graph element, we have to supply three parameters to the S-function. They are: • The bond graph element name.
• Whether the element is a decider bond or not for the producer -consumer correspondence to the input -output port to be decided. • The name of the junction to which it is connected.

S-function library for bond graph elements
Bond graph modeling is an emerging field especially in electrical, mechatronics and electro-mechanics, and a bond graph engineer may feel the necessity to define his own element while making a model. With the set of callback routine and function available to MEX files, any complex constitutive equation can be written for a new element. To provide support for complex variables and vectors in Bond Graph, a tool library using C-MEX S-functions with data propagation through shared memory, is developed. The MEX file for the library has been written in C/C++. After compiling and debugging, C/C++ MEX S-Function are masked with bond graph icons to distinguish between different elements.
Each element of the bond graph library has two common input/output blocks along with with a middle block (Fig.11). The code for the middle block is specific to the element it implements. After placing the three S-functions block in the subsystem, the subsystem is masked. The element's mask screen has a help at the top and parameter entry text boxes, below. There is a check box to specify whether the bond is decider (Fig.12). The parameters entered in the mask screen are manipulated by the S-functions underneath to initialise the element before the simulation cycle starts.
The capability of S-functions to support arbitrary input dimensions is exploited in the tool box. The actual input dimensions can be determined dynamically when a simulation is started by evaluating the dimensions of the input vector driving the S-function. This feature allows the same element to handle a scalar or a vector input as the case may be, without declaring it apriori. The library is available in the standard format of simulink (Fig:13). The required elements can be had for Pick and Place from library after navigating down (Fig:13(a) and Fig:13(b)). The tools available in the mask's -icon graphics support, is utilised to give a natural iconic representation to the element subsystem.

Examples of tool box
Using the tool box, the circuit in Fig.14(a) is modeled in bond graph ( Fig.14(b)). The simulation results are given in Fig.14(c) and Fig.14(d). As can be observed the library is able to handle complex quantities quite accurately. The circuit is modeled in bond graph using the switched junction as shown in Fig.14(e). The simulation result of the effected state variable is shown in Fig.14(g).

Simulation results for IM model
A rotating electrical machine can be viewed as a machinery which converts one form of energy into another. More specifically it converts electrical energy into mechanical energy or vice-versa. Magnetic energy is used as a conversion medium between electrical and mechanical energy.

Axis rotator element
This generalised concept for electrical machine modeling needs 'Axis Rotator', a new bond graph element (Umesh Rai & Umanand, 2008;2009a; to mathematically model a electrical commutator (Fig.15). The constitutive relationship for the flow and effort in the bonds is given by Eqn.
(3).  cos α (i, k) refers to the spatial angle between the k th winding's axis and the i th winding axis with respect to the bond under consideration. Λ m is the mutual permeance (inverse of reluctance) of the magnetic core, P m is the reactive power required to magnetise the core.

Induction motor model
A bond graph model of 3φ doubly fed induction motor using the Axis Rotator element is shown in Fig.16. There are six sets of electric energy input ports, three each for stator and rotor, in the model. The motor shaft represents a mechanical output port. The air gap is represented by the AR with six connection bonds terminating at it, each representing a set of    winding. The permeance parameter of the AR represents the mutual coupling effect of the flux. The iron loss is represented by the dissipative port. The stator and rotor energy ports are described with different set of port parameters.
The dissipative, entropy producing fields, R as , R bs , R cs , R ar , R br , R cr and R i are non-equal non-linear resistances depending on the underlying physical system. R i represents the iron losses of the core. In a similar manner, the model represents the permeances L asl , L bsl , L csl , L arl , L brl , L crl and L m . There is no linearity restriction on the above parameters. They could be constants, functions or even a lookup table, without loss of generality. The value of the lumped parameters are different for different phases, dependent on the electric and the magnetic energy they represent. Similarly the energy sources feeding the different windings are represented by s 1 , s 2 and s 3 . For balanced supply voltages the voltage peaks and frequency would be same, with a phase difference of (2π/3) to one another. The three domains of electrical, magnetic and mechanical are clearly brought out in Fig.16.
At the shaft the developed electromagnetic torque as a function of the stator, rotor currents and the angle between them is represented as an effort source. The electromagnetic torque provides the effort at the one junction against the inertial, frictional and load torque components. The feedback information on flow at this junction, which gives the measure of rotor speed is transmitted to the AR for calculating the instantaneous angle of the rotor windings.

S-function model of induction motor
The power of S-function is demonstrated by firstly implementing the complex AR element using C Mex S-function and making it a part of bond graph library. As discussed in the above section, there is no linearity or balance supply constraint on the model. The increased complexity can easily be handled by the bond graph library. The causal model implementation of squirrel cage induction motor is shown in Simulink (Fig.17). The machine starts from stall. A step load is applied to the motor at 0.5sec. The simulation results of speed curve for various step load are presented in Fig.18(a). Similarly the current curves and the torque curve for a specific load are at Fig.18(c) and Fig.18(e) respectively. The transition of rotor currents to slip frequency can be distinctly seen at Fig.18(c). The simulation results of the bond graph model in S-function co-related well with the speed curves obtained by d-q model of the induction motor implemented by functional block in Simulink ( Fig.18(b), Fig.18(d) and Fig.18(f)).

Conclusion
The power of S-function in customising MATLAB/Simulink ® environment suitable for a specific modeling need is illustrated in this chapter. Two different approaches for implementing a customised Simulink element is then discussed with their advantages and disadvantages. Later the Bond graph approach of modeling is briefly introduced. Level 2 C MEX S-function technique is then used to develop a library of bond graph element.
This library of bond graph elements can handle both the scalar and vector or complex variables without declaring apriori, a distinct advantage. A new element -Axis Rotator, used for representing rotating magnetic field is included in the library. The ability to handle complex variables along with the rotation enables the elements in the library to be used for modeling rotating frames as that existing in an electric machine. The library also incorporates switched junctions, which allows for the modeling of switches in any circuit. In developing the library, the Shared Memory Concept is used. By using shared memory concept, the memory requirement comes down as only the pointer to the memory location are passed and not the data values.