Component-oriented acausal modeling of the dynamical systems in Python language on the example of the model of the sucker rod string

Typically, component-oriented acausal hybrid modeling of complex dynamic systems is implemented by specialized modeling languages. A well-known example is the Modelica language. The specialized nature, complexity of implementation and learning of such languages somewhat limits their development and wide use by developers who know only general-purpose languages. The paper suggests the principle of developing simple to understand and modify Modelica-like system based on the general-purpose programming language Python. The principle consists in: (1) Python classes are used to describe components and their systems, (2) declarative symbolic tools SymPy are used to describe components behavior by difference or differential equations, (3) the solution procedure uses a function initially created using the SymPy lambdify function and computes unknown values in the current step using known values from the previous step, (4) Python imperative constructs are used for simple events handling, (5) external solvers of differential-algebraic equations can optionally be applied via the Assimulo interface, (6) SymPy package allows to arbitrarily manipulate model equations, generate code and solve some equations symbolically. The basic set of mechanical components (1D translational “mass”, “spring-damper” and “force”) is developed. The models of a sucker rods string are developed and simulated using these components. The comparison of results of the sucker rod string simulations with practical dynamometer cards and Modelica results verify the adequacy of the models. The proposed approach simplifies the understanding of the system, its modification and improvement, adaptation for other purposes, makes it available to a much larger community, simplifies integration into third-party software.


INTRODUCTION
As is known, component-oriented simulation modeling is based on the separation of a complex system model into simple components. The component describes the mathematical model of the corresponding physical object (mass, spring, electrical commercial simulation environments for Modelica-OpenModelica, JModelica.org, Wolfram SystemModeler, SimulationX, MapleSim, Dymola, LMS Imagine.Lab AMESim. The known problem of specialized languages is the complexity of the modifications and improvements and a relatively small community of developers. They are not very well suited for experimenting with evolutions of modeling capabilities (Elmqvist, Henningsson & Otter, 2016). In particular, developers of some languages have encountered the problem of variable structure systems modeling where the structure and number of equations can change at run-time (Fritzson, Broman & Cellier, 2009;Nikoli c, 2016). Some problems can be solved by using interfaces to general-purpose languages (Åkesson et al., 2010;Hedengren et al., 2014). But it is usually more difficult to learn a new language than to learn a component or a library of a familiar programming language.
As a rule, general-purpose languages, in comparison with specialized languages, are more widespread, easy to learn thanks to typical imperative and object-oriented constructs, have wider applicability, better interoperability with the third party software and a large number of heterogeneous packages. Therefore, the mentioned problems are less common in modeling systems that are based on general-purpose programming languages: PyDSTool (Clewley et al., 2007)-Python-based Dynamical Systems Toolkit with support for symbolic manipulation, hierarchical structures and hybrid models; Ariadne-a C++ library for formal verification of cyber-physical systems, using reachability analysis for nonlinear hybrid automata (Benvenuti et al., 2014); Assimulo-Python-package that combines a variety of different ODE/DAE solvers via a common high-level interface (Andersson, Führer & Åkesson, 2015); DAE Tools-equation-based object-oriented modeling, simulation and optimization software with hybrid approach (Nikoli c, 2016); Modia.jl-Modelica-like language that is directly defined and implemented with Julia's meta-programming constructs and is designed tightly together with the symbolic and numeric algorithms (Elmqvist, Henningsson & Otter, 2016); SimuPy-a Python framework for simulating interconnected dynamical system models (Margolis, 2017); Sims.jl-a Julia package for equation-based hybrid modeling and simulations, similar to Modelyze (Short, 2017); GEKKO-a Python package for machine learning and optimization of mixed-integer and differential algebraic equations (Beal et al., 2018). However, most of these systems either have a complex code that is difficult to understand and modify (e.g., have their own symbolic processors), or use state-of-the-art ODE/DAE solvers that are implemented in low-level languages, which rarely allow modification to an untrained user. The implementation, modification and improvement of such systems can be simplified if the difference equations are used to describe the model instead of differential equations. In addition, difference equations are also often used to model dynamical systems. In general, the modeling system should allow various types of equations.
The advantages of modeling systems based on general-purpose programming languages are described in detail in papers (Nikoli c, 2016;Elmqvist, Henningsson & Otter, 2016). Python language (Van Rossum & Drake, 1995) is a good choice mainly due to its features: multi-paradigm, object-oriented, intuitive with code readability and improved programmer's productivity, highly extensible, portable, open-source, large community and extensive libraries as mathematical libraries SymPy and SciPy. SymPy is a Python library for symbolic mathematics (Meurer et al., 2017). SciPy is a fundamental library for scientific computing (Oliphant, 2007).
This work suggests the principle of developing simple to understand and modify Modelica-like system based on the general-purpose programming language Python. The principle consists in: (1) Python classes are used to describe components and their systems, (2) declarative symbolic tools SymPy are used to describe components behavior by difference or differential equations, (3) the solution procedure uses a function initially created using SymPy lambdify and computes unknown values in the current step using the known values from the previous step, (4) Python imperative constructs are used for simple events handling, (5) external DAEs solvers can optionally be applied via the Assimulo interface, (6) SymPy package allows to arbitrarily manipulate model equations, generate code and solve some equations symbolically. The principle of the system is described by examples of models of a sucker rod string that is used in the oil industry to connect surface and downhole components of a rod pumping system.

Description of modeling principles in Modelica language
First, we describe the modeling principles in Modelica using an example of a simple mechanical translational oscillator. The oscillator consists of such components as Mass, SpringDamper and Fixed (Fig. 1). The SpringDamper component is designed to simulate the elastic-damper properties of the damped oscillator. The Mass component simulates the inertial properties of the oscillator. The Fixed component simulates the fixed point of the oscillator. The module code that describes this model is explained below (Listing S1). In order to simplify the model, these classes differ slightly from the corresponding classes of the standard Modelica library (Fritzson, 2015).
A class in Modelica describes the set of similar objects (components). The Flange class describes the concept of a mechanical flange. Its real-type s variable corresponds to the absolute position of the flange. Its value should be equal to the value of the s variables of other flanges connected to this flange. The real-type f variable corresponds to the force on the flange. It is marked by the flow keyword, which means that the sum of all forces at the connection point is zero. The Fixed class describes the concept of a fixed component with one flange, for example, fixed1 (Fig. 1). It has the real-type s0 variable, which corresponds to the absolute position of the flange, and the flange object of the Flange class, designed to connect this component to others. The s0 variable is marked by the parameter keyword, which means that it can be changed only at the start of the simulation. After the equation keyword, an equation describing the behavior of this component is declared-the flange object position must be equal to the s0 value. The Mass class inherits the Transl class and describes the sliding mass with inertia. The example of such component is mass1 (Fig. 1). The extends Transl command means inheriting members of the Transl class in such a way that they become members of the Mass class. That is, the Mass component will also have two flanges (flange_a and flange_b). In addition, this class has the m parameter (mass) and the variables: s (position), v (speed), a (acceleration). The expression start=0 is the default initial condition. After the equation keyword the system of the differential and algebraic equations, which describes the behavior of this component, is given. The der keyword means the derivative with respect to time t (v = ds/dt, a = dv/dt). The SpringDamper class inherits the Transl class and describes the linear 1D translational spring and damper in parallel (Listing S1). The example of such component is springDamper1 (Fig. 1). The class has the parameters c (spring constant), d (damping constant) and the variables s_rel (relative position), v_rel (relative speed), f (force at flange_b). After the equation keyword the system of differential-algebraic equations of this component is given.
The Oscillator class describes the spring-mass system (Fig. 1). It contains three components mass1, spring1, fixed1, which are described by the classes Mass, SpringDamper and Fixed, respectively. The values of parameters and initial conditions of these components are shown in round brackets. The model code can be prepared using any text editor or the Modelica Development Tooling module (Pop et al., 2006) of the Eclipse development environment. Simulation of a model requires the OpenModelica environment (Fritzson et al., 2005). To start calculations and to plot the curve, which describes the position of mass1 component with time, enter the following into the OpenModelica Shell: loadModel(Modelica) loadFile("Pycodyn.mo") simulate(Pycodyn.Oscillator, stopTime=10) plot(mass1.s) The typical stages of translating and executing a Modelica model are (Fritzson, 2015): translation (obtaining a flat set of equations, constants, variables and function definitions from Modelica-code), analysis (equations sorting, convert the coefficient matrix into block lower triangular form), optimization (elimination of most equations, converting equations to assignment statements), code generation (obtaining a C-code), compilation (obtaining an executable) and simulation. However, there are alternative ways of executing Modelica (Fritzson, 2015).

Description of modeling principle in Python
The principle of component-oriented modeling in Python is described below and an example of the implementation of a similar oscillator model is shown.
1. Components are described by Python classes that are structurally similar to Modelica classes and have the following attributes: constant parameters and SymPy symbols (analogs of parameters and variables in Modelica), SymPy symbolic equations (difference or DAE), pins for connecting components into a system (analogs of connectors in Modelica). A dynamic system consisting of components is also described by a Python class which attributes are a list of components and equations (together with additional equations for connecting components).
2. If the dynamic behavior of the components is described by difference equations, then the user must describe these equations in the class by replacing the derivatives with the selected difference scheme (e.g., by the Euler method).
3. The initial conditions are substituted into these difference equations and unknowns are found by solving a system of nonlinear equations at each step or using a function that was initially created using the SymPy lambdify function (translates a SymPy expression into an equivalent numeric function) and calculate unknown values at current step from the known values from the previous step without the need to solve the equations.
4. At each step, the if statement checks for discrete events that depend on state or time. During event handling, initial conditions, components, or equations can be changed.
5. If the dynamic behavior of the components is described by DAEs, then the Assimulo interface to the DAEs solvers is used, which has an effective discontinuity handling procedure.
6. The SymPy package allows arbitrary manipulation of model equations and code generation. You can solve some algebraic or differential equations symbolically. The DAE system of equations should be simplified, transformed into an ODE, and solved by the SymPy dsolve function. Now we will develop the pycodyn module with similar components in Python (Listing S2). The behavior of the components will be described using the difference equations. For simplicity, we will use the Euler method. As a result, the system of components connected by flanges will be described by the system of the difference equations.
First, import the sympy module and the standard mathematical module math. It is important to distinguish the functions of these modules. from sympy import Ã import math Create the global variable dt (time step).

dt=0.1
If you only need to obtain the system of equations in a symbolic form, then this variable must be an instance of the Symbol class from the sympy module: Translational1D is the basic class of mechanical 1D components with translational motion. The __init__ method is called when an object of this class is created and has two parameters-the name of the component (name) and the dictionary of its attributes (args). For component attribute naming, we use the following notation. At the beginning of the name, the x, v, a, f symbols mean position, speed, acceleration and force, respectively. At the end of the name, the p symbol means the value at time t-dt.
The Mass class describes the mass concentrated at a point, which has a translational motion. It inherits Translational1D class. The __init__ constructor calls the constructor of the base class Translational1D and sends to it the parameters name and locals(). The latter is a dictionary of local variables self, name, m, x, xp, v, vp, a, f1, f2.
class Mass (Translational1D): def __init__(self, name, m=1.0, x=None, xp=None, v=None, vp=None, a=None, f1=None, f2=None): Translational1D.__init__(self, name, locals()) # base class constructor call The behavior of this component is described by a system of equations self.eqs. For example, for an m1 component: A list of additional equations can be generated for each component flange using the pinEqs method described above. The first element of the self.pins list is the dictionary dict(x=self.x, xp=self.xp, f=self.f1), which means that the x, xp positions on the flange will be equal to the self.x, self.xp attributes of this component, respectively, and the force f on the flange will be equal to the self.f1 attribute. The same applies to the second element of the list.
The SpringDamper class (Listing S2) describes the translational 1D spring and damper, which are connected in parallel. It inherits Translational1D class. In addition to the attributes described above, it has the following attributes: spring constant c, damping constant d, relative velocity between flanges vrel. The behavior of this component is described by a system of equations self.eqs. E.g. for an s1 component: This component also has two flanges and it is possible to generate a list of additional equations using the pinEqs method.
The Force class (Listing S2) describes a 1D force with a translational motion of the application point. The value of the f force can be constant or variable. It inherits the Translational1D class and has one flange.
The System class (Listing S2) describes the system of components connected by flanges. The constructor __init__ gets two parameters-the list of components els and the list of additional equations eqs, which usually are created using pinEqs method. The system components are stored in the self.els list and the self.elsd dictionary. The list self.eqs contains all system equations and is created by joining the equations of all components with additional equations eqs.
The solve method of this class solves a stationary problem. It returns the solution of a system of equations with conditions ics-a dictionary with known values of variables. To solve a system of equations, it can use the SymPy solve function, but its algorithm is very slow. It is possible to use fast algorithms for solving equations, for example, the function scipy.optimize.root from the SciPy library, which supports many effective methods for solving nonlinear systems of equations. In this case, the call of the SymPy function solve(eqs) must be replaced with the call of the self.solveN(eqs)method, which adapts the system of equations for SciPy and solves it using scipy.optimize. root.
The solveDyn method solves a non-stationary problem. It receives three parametersthe dictionary with initial state state, the final time value timeEnd and the fnBC function that returns the dictionary to update the state. First, the time variable t is assigned an initial value. In the while loop with the condition t<timeEnd, the following instructions are executed: previous step variables (xp, x1p, vp, etc.) are assigned the values of the initial state state, the values of the boundary conditions are updated, the system of equations is solved by calling the self.solve method, solutions are assigned to the dictionary state, the results are saved, the time value increases by dt. After the loop is completed, the method returns the results as T and Res lists. These results can be represented in the form of plots using the matplotlib library. But the use of the self.solve method can be acceptable only for very frequent discontinuities that require the re-creation of equations. Since at each step this method creates and solves a system of nonlinear equations, the calculations can be very time-consuming. In most cases, it should be replaced by the solvN method, which at each step finds unknown values by passing the values found in the previous step to the ceqsf method. At the beginning of the simulation and after the discontinuities, this method must be created using SymPy lambdify function that transforms SymPy expressions to lambda functions which can be used to calculate numerical values very fast. This is done in the createCurEqs method.
Event processing is performed at the end of each step by calling the user-defined event handling method event(state). In it, the if statement checks a determined condition with state. If the result is True, then the event is handled, for example, new boundary conditions are created and createCurEqs is called. You can easily implement modeling of variable structure systems by calling in the event method the constructor of the System class (with new values of els, eqs) and the createCurEqs method.
If differential equations are used in the components, then the functions and their derivatives should be distinguished. The presence of the symbol "D" in the name of the variable means derivative. For example, m1_Dx is a derivative of m1_x. Class descriptions of such components will be more like Modelica classes (Table 1). In this case, to solve the DAEs in the form 0 = F(t, y, y'), the solveDAE method from the pycodynDAE (Listing S3) module is used. Assimulo was used as an interface with ODE/DAE solvers such as SUNDIALS IDA (Hindmarsh et al., 2005) or DASSL (Petzold, 1982). The solveDAE method forms the residual method and initial values for the time, states and state derivatives required by a DAEs solver. The residual method takes as input time t, state y, state derivative y' and returns a residual vector (zero if a solution is found). Argument lists for the residual are prepared by the residualArgs method. It is also possible to create user-defined functions state_events and handle_event for event tracking and handling in discontinuous problems for Assimulo. Some problems in pycodyn, which is formulated using difference or differential equations, can be solved symbolically using the SymPy functions solve and dsolve. In particular, the solve function, which symbolically solves equations and systems of equations, helps to form the ceqsf method mentioned above. The dsolve function solves any supported kind of ODEs. Therefore, DAEs needs to be transformed into ODEs by simplification.

USE CASES
For testing purposes of pycodyn, we consider models of sucker rod strings. Let's take a look at the steel sucker rod string, in which the length is 1,510 m. Such a string and its practical dynamometer card are described in (Belov, 1960). The upper section of the string consists of 695 m rods with a diameter of 22 mm, and the lower section consists of 815 m rods with a diameter of 19 mm. This string will have a total mass of 3,961 kg, a total weight in the liquid of 34,687 N, a spring constant of 44,650 N/m, a damping constant of 2,121 N•s/m. Liquid weight above the pump with a diameter of 43 mm will be 18,499 N.

Use case 1: simulation of free vibrations of the sucker rod string
First, simulate the free vibrations of the string using the Modelica language. Initially, the string is stretched by moving the lower end of the string (mass1.s) by one m. After the lower end is released, free vibrations will begin. Use the model (Listing S1) with parameter values: mass1.m=3961.0, spring1.c=44650.0, spring1.d=2120.7, and with initial conditions: mass1.s=-1, mass1.v=0. Simulation options: stopTime = 10.0, numberOfIntervals = 500, tolerance = 1e-006, method = 'dassl'. Now perform the simulation of free vibrations of the sucker rod string in pycodyn (Fig. 1). In the separate module (Listing S4) create the components: spring-damper s1 and mass m1. In round brackets, there are the values of the attributes-the name and the known parameters.
To build the single-section model Pumping in Modelica (Listing S1), use the components of the oscillator and additional components: motion1 (describes the movement of the upper point) and force1 (describes the forces acting on the lower point). This is a simplified model that does not take into account other types of loads (Kopey, Kopey & Kuzmin, 2017). The stroke length of the upper point is 2.1 m, the number of double strokes per minute is 6.4. During the downstroke, the weight of the rods acts on the lower point. During the upstroke, the weight of the liquid is added to it. This is shown in the algorithm section of the Pumping model: The analog of this single-section model in Python is shown in (Listing S8). Now in the new module (Listing S9) create the two-section Python-model of the sucker rod string, in which each section of the string is modeled by three components. The model of each section consists of three 1D mechanical translational components: SpringDamper, Mass and Force (Fig. 2)  Create the components: the spring-damper of the first (upper) section s1, the mass of the first section m1, the weight of the first section f1, the spring-damper of the second section s2, the mass of the second section m2, the weight of the second section with the weight of the liquid f2.  The complete list of equations for this system s.eqs in the SymPy format: Solve the static problem-the string under the maximum static loads. A=2.1/2 # amplitude n=6.4/60 # frequency return A Ã math.sin(2 Ã math.pi Ã n Ã t) # position The force function returns the value of the force on the pump plunger F, depending on the value of its speed v. If the speed is less than zero (downstroke of the string), the function returns the weight value of the second section. Otherwise, the function returns the sum of the second section weight and the liquid weight above the plunger. This function should be smoothed when the sign of the velocity changes, for example, using the math. tanh hyperbolic tangent function. The simulation of the variable structure system (breakage of the second section) by the Euler method (dt=0.1) is implemented in (Listing S10). If the force at the top of the second section is greater than 56,000 N, then the section breaks off. The user method event is created to handle the event. At each step, this method checks the condition state[s1.f1]>56000. If the result is True, then an event occurs. The handling of this event consists in changing the components of the system (only s1, m1, f1 remain after the breakage), changing additional equations at the connection points, changing the boundary conditions fnBC (the weight of the second section and the weight of the liquid are zero). The ceqsf method is also updated by calling createCurEqs. Analogous single-section (Listing S11) and two-section (Listing S12) models, based on differential equations, are developed using the pycodynDAE module. The DAE system should be supplemented with equations that describe the position of the upper point and the force that acts on the lower point. For example: eq=s.eqs+Tuple(Eq(s1.x1, A Ã sin(2 Ã pi Ã n Ã t)), Eq(m1.f2, Piecewise((fs, m1.v<0), (fs+fr Ã tanh(abs(m1.v)/0.01), m1.v>=0)))) Alternatively, you can replace the symbols s1.x1 and m1.f2 with the right-hand sides of these equations.

SIMULATION RESULTS
The results of the simulation of free vibrations are shown in Fig. 3. The frequency of free vibrations corresponds to the theoretical natural frequency of the harmonic oscillator v = (j/m) 1/2 = (44,650/3,961) 1/2 = 3.357 rad/s, where j is the stiffness, m is the mass. Such vibrations occur during normal operation of the pump due to the sharp removal or application of the load (the pump valve opens or closes) and are noticeable in the upper and lower parts of the dynamometer card (Belov, 1960). The differences are explained by the use of unequal difference schemes. The results obtained analytically and by IDA/DASSL solvers are almost equal, therefore they are shown by a single curve in the figure. The global error values (at t=1 s) for the Euler (dt=0.1), trapezoidal rule (dt=0.1), Euler (dt=0.01), DASSL (OpenModelica), SUNDIALS IDA methods are, respectively, 0.344, 0.034, 0.046, 6.35E-05, 1.46E-05. The simulation time is 0.25, 0.5, 2.46, 0.29, 0.028 s, respectively. The total time of the module execution (total simulation time for OpenModelica) is 4.1, 5.5, 6.3, 5.1, 1.9 s, respectively. These data were obtained for this configuration: simulation interval 0-10 s, CPU 2.5 GHz, Python 3.7, NumPy 1.16.4, OpenModelica 1.12, Sundials 2.6.
The results of the pumping process simulation (Fig. 4) correspond to practical dynamometer cards obtained on real wells (Belov, 1960). Single-section (Listing S8) and two-section (Listing S9) models, simulated by the Euler method (dt = 0.1), somewhat distort the right and left sides of the card and smooth the upper and lower sides (Fig. 4A). Single-section (Listing S11) and two-section (Listing S12) models, simulated by the IDA, give somewhat larger values of maximum load and lower values of minimum load (Fig. 4B). The two-section model is more adequate.
The dynamometer card of the string breakage model (Fig. 5) corresponds to practical dynamometer cards with their typical flat shapes (Belov, 1960).

DISCUSSION
It is noticeable that the global error for the Euler method (dt = 0.1) is much larger, therefore it is not advisable to use it for the free vibration problem. In addition, numerical solution procedure in pycodyn has a low performance if difference equations are used. In the future, it is planned to improve performance, for example by using Cython.
The differences with the practical dynamometer card are explained by the fact that the real string has a greater number of degrees of freedom. For a more adequate simulation, you need to increase the number of sections of the model (Kopey, Kopey & Kuzmin, 2017) or use the wave equation, which is a partial differential equation (Gibbs, 2012). In addition, it is difficult to know the exact value of the damping constant, which depends on many factors (Kopey, Kopey & Kuzmin, 2017). In general, all these models can be used for approximate modeling of the pumping process.
The Python language allows a simple modification and improvement of pycodyn. In the future, it is planned to extend the set of the components (e.g., create electrical and hydraulic components), develop support for hierarchical subsystems and the tools for building models using component diagrams. To implement hierarchical subsystems, you can move functions for solving equations to a separate module and add pins to the System class. The problem in the form of difference equations is usually more difficult to formulate. But SymPy can be used to automate the conversion of differential equations to difference equations. In order to simplify the code of the pycodyn and pycodynDAE modules, the algorithms for equations sorting, eliminating, and simplifying are not implemented. This is planned to be implemented in the future using SymPy.

CONCLUSIONS
The Python-classes that allow creating the Modelica-like models in Python without the need to study and apply specialized modeling languages are developed. The suggested approach simplifies the understanding of the system, its modification and improvement, adaptation for other purposes, makes it available to a much larger community, simplifies integration into third-party software. Difference or differential equations can be used to describe components. Using difference equations and the pycodyn module allows simplifying the implementation of the hybrid modeling, variable structure systems modeling and the requirements for the modules for symbolic mathematics and for solving equations. It is also well suited for experimenting with evolutions of modeling capabilities. In particular, one can experiment with arbitrary difference schemes, make arbitrary symbolic manipulations and modify numerical solution procedures. However, the pycodynDAE module can provide higher accuracy and performance using third-party DAEs solvers that are suitable for stiff problems. With SymPy, some tasks can be solved symbolically.
The comparison of simulation results of sucker rod string with practical dynamometer cards and Modelica models verify the adequacy of the models. The pycodyn framework can be used to study the principles of component-oriented modeling and for various kinds of experiments on its new features. The source code is freely available under the GNU GPLv3 open-source license from the GitHub (https://github.com/vkopey/pycodyn).

ADDITIONAL INFORMATION AND DECLARATIONS
Funding