DynSim

DynSim2, the beginnings of a second version of DynSim, is a simple process control (and dynamic systems) simulator in Python. DynSim, the original, was written in C and designed for speed, but was not scriptable and lacked some of the flexibility of a version written in Python. To some extent, the speed of modern computers can make up for the loss in speed from using an interpreted language. However, any hard number crunching elements (spectral analysis, NNs etc.) will still use the compiled extensions to Python.

Unfortunately, DynSim now appears to be a registered trademark of Invensys, so a name-change for this particular DynSim may be in the winds.

Features

Very lightweight and fast simulation classes.
A class library concentrating on practical process control algorithms:
- control (basic modulating control building blocks which work similarly to those found in PLCs and DCSs)
- process (control oriented process classes, including representation of 2-phase flow streams)
- mechanics (basic mechanical systems - barely started)

Motivations and Philosophy

For experienced model developers, a good code-based simulator can be (in the opinion of the author):
- more concise,
- more readable / understandable, and
- more powerful, than a graphically constructed simulation tool (such as Simulink), especially if dealing with non-trivial one-off models.
Provide the simplest possible continuous- and discrete- time simulation language that will encourage the use of simulation in the most trivial of problems (rather than the hurdle causing simulation to be reserved for only the most complex of situations).
Exploit one of the most friendly object-oriented languages to make simulation construction the easiest possible, but have a language which can be taken all the way to include phenomenological models, nice GUIs, databases and web interaction as desired.
Choose a language which is noted for its ability to integrate third-party models and algorithms (be they fortran, C or C++ or others) by being able to quickly write wrappers for these modules.
Don't create a new, simulation-specific language. Just extend an existing language with object-oriented elements (object classes, including operator-overloading methods).
Provide the essential building blocks (classes) from which to create the simulation, with an emphasis on conventional and advanced process control functionality.
The blocks should always have sufficient initialisation functionality, to allow initialisation (e.g. manual-override) during a simulation, as would be done with "house-keeping" logic in a real process control scheme implementation. This tool is for real-life control design, not simplified academic tutorial.
Within the same language environment, have the full Scipy library of advanced matrix and numerical methods to include in the simulation and wrap around the simulation as necessary (e.g. solving a simulation for best-fit parameters).
Provide a set of simple examples to guide the first-time users.

Graphical vs Code

Conciseness, readability and power are all obviously in the eyes of the beholder. However, for object oriented systems, where reuseable blocks/classes/objects are glued together with the "custom" elements of your simulation, it is common for these glue elements to be unconventional, complex or special purpose. Often they represent complex equations or algorithms. Such equations and algorithms are almost always easier to program, and more readable, in code form. Obviously a good object oriented language like Python allows even more reuseability through inheritance by being able to take an existing class and specialise it to suit your purpose. Have you ever used a graphical language (with X, +, - blocks etc.) to program a mathematical equation? What might start out as an obviously recognisable equation form in your expertise domain, turns into a rabbit warren of blocks and lines. In Simulink, it seems better to resort to embedding some m-file code into your simulation when you need to put in an equation (even accounting for the performance hit).

Graphical model construction is somewhat misleading to casual users/observers in that it appears that you can construct complex models without programming any code. Of course, this is rubbish! It is just a different form of code; one that still needs to be learnt in great detail (like a normal programming language). However, in order to pack the programming into this (perhaps "unnatural") graphical form, much of the flexibility of a normal programming language is lost. This usually requires the further learning of tricks and techniques to get around the built-in limitations of the graphical construction.

Integration Method

It might be argued that you need "proper" Runge-Kutta etc integrators (which are somewhat problematic for block- or object- oriented modelling) to get good fidelity from simulations. However, more often than not, Process Control models and simulations are often taken overboard with accuracy and detail, with combined controller and process differential equations being solved quite accurately. However, if your objective is to simulate more realistic scenarios, then as an example, the discrete-time implementation of PID controllers (in PLCs and DCSs) needs to be simulated with Euler integration with fixed time-steps of the order of 0.5 seconds or 1.0 second. In this case, differential equation equivalents are inaccurate and the simulation has lost most of its motivation to provide hi-fidelity solutions to the differential equations. This philosophy is exploited in DynSim2 with only manually-set, fixed time-steps being used through Euler or trapezoid integration (built into each simulation class); This helps to "keep it simple if possible, but make it possible to do difficult things" (e.g. through fine timesteps or multi-rate timesteps).

Another advantage of fixed-timestep integration is that the set of calculations done during a simulation is fixed (apart from any sequencing elements). This is advantageous (and almost a pre-requisite) for enclosing the simulation in an optimisation loop, such as would be the case for fitting model parameters to plant data and optimising tuning parameters for controllers. In these cases, variable-step solvers provide an inconsistent result change for a very small change in parameters, due to the number of time-steps for the simulation changing - essentially it generates noise in the objective function of the optimiser. Optimisers relying on smooth objective functions will not work well with this noise.

Getting Started

Go to the Tutorial page for a complete walk-through of constructing, executing and analysing a simulation.

Documentation

The Documentation page provides the full documentation of DynSim, its usage, the class library, and advanced topics.

A Simple Example

A couple of these classes can be quickly put together into a simulation of a PID control loop with a process represented by a first-order-plus-deadtime (FOPDT) characteristic using the following example script (from the source file):

"""Simple example of dynsim use for a control loop.
Simulates and plots control of a first-order plus deadtime (FOPDT) process
with a PID controller.  100 seconds at stepsize=0.1 second are simulated."""

from pctools.dynsim import control

# create the simulation objects or "blocks"
ctl = control.PID(P=1.8, I=5.0, D=0.0)  # make a PID controller
plant = control.FOPDT(gain=1.0, delay=2.0, tc=5.0) # first-order plus deadtime model
record = control.Log(name='Logged Data', legend=('PV', 'CO', 'SP'))  # set up a logger


PV = 0.0  # provide PV for first time-step (PV not yet generated by plant)
SP = 0.0   # initial value of Setpoint before step test

# simulate
for t,dt in control.timesteps(finish=100.0, timestep=0.1, plot=True):
    if t>10: SP = 40.0  # step change in setpoint
    out = ctl(SP, PV)  # run the PID controller
    PV = plant(out)  # run the process dynamics
    record([PV, out, SP])  # record a log of the main variables

This script runs 1000 time-steps in about 300 milliseconds and produces the following plot window at the end.

If you're thinking "Gee, that seems like a lot of code to do that!". However, if you look carefully, more than half of the code-text is comments. Here is the code again without the comments - its much cleaner and concise, but adding comments still wins for good understanding of the code.

from pctools.dynsim import control

ctl = control.PID(P=1.8, I=5.0, D=0.0)
plant = control.FOPDT(gain=1.0, delay=2.0, tc=5.0)
record = control.Log(name='Logged Data', legend=('PV', 'CO', 'SP'))

PV = 0.0
SP = 0.0
for t,dt in control.timesteps(finish=100.0, timestep=0.1, plot=True):
    if t>10: SP = 40.0
    out = ctl(SP, PV)
    PV = plant(out)
    record([PV, out, SP])

Download it

You can find the latest version of DynSim in the Weblog below. Look for the latest entry with the title "Uploaded new version ...". The zip file will contain all of files you need. Find your installed version of Python and navigate to "...\python\Lib\site-packages". Unzip the dynsim directory there. That's it! You should now be able to start up Python and from dynsim import control etc.

Application Ideas

Simulation of control schemes plus process dynamics (the main idea)
PID controller tuning optimisation (for a known FIR or FOPDT response curve)
fitting of a FOPDT or FIR model to logged response data
- based on both CO and PV logs
- allow fitting of a PRBS (or arbitrary excitation) response
fitting of a FOPDT or FIR to closed-loop response (assumes PID form is correct for at least the terms involved)
NB: These may require model to be constructed as a class (for embedding into an optimisation loop).

ToDo

more class call parameter checking (type, range) to help avoid exceptions inside dynsim code (hard for users to debug cause)
- use exceptions with object "name" attribute, parameter, and fault to inform the user.
- or perhaps use "assert" with the appropriate message (will this generate a call stack to track down where, or don't we want this since it will lead into dynsim code?)
Add a PID loop example with post analysis of data to show meaning of autocorrelation, crosscorrelation, spectrum, system-ID, ...)
- test bed for pctools and interpretation tutorial
Add demo of tuning a PID loop with MatPlotLib widgets (can use widgets with pylab!)
Decide on a convenient method to add GUIs to real-time plots and/or GUI around repetitive simulation (such as a PID tuning simulator).
- Qt-Designer option
  - allows complete customisation of GUI using the GPL designer to produce a .ui file which can be imported into a python script
  - create a class which accepts a dictionary of gui_item_name to variable relationships that automatically hooks up the GUI components to the variables in a simulation.
  - main issue is event loop issues between real-time plotting and real-time GUI
- add a widgets class to dynsim.control
  - separate window to start with, but option integaration with log through log.fig
  - configure with a list or dictionary of widgets
  - automatically add space and scale the widgets
- Add a module which allows Vpython displays (for real-time interactive simulations). Its already used for the CartPole demo, but this is a window specific to the simulation object rather than a GUI display. Try to expand it to become a general purpose GUI for the whole simulation.
  - for animated outputs, but especially ...
  - for (button, slider, menu) controls in a window
- or, use MatPlotLib's controls for real-time input
  - however, this seems difficult to integrate with a running simulation because its event loop is in show(), whereas Vpython's is external and can be integrated.
  - MatPlotLib's interface would be good though, since the other pylab functions are already used and the controls can be placed around normal plot axes, and it can automatically be integrated with pydee.
  - http://www.scipy.org/Cookbook/EmbeddingInTraitsGUI has an example of Matplotlib using Enthought's Traits with wx.
Consider encapsulating models in a "model" class (either as standard, or for special purposes):
- to enable a GUI simulation control (is this useful, given that script cannot be changed between GUI "runs"?)
- to enable "automatic" optimisation functions (to fit parameters, optimiser PID tuning, etc)
- "model.step" method would only do one time-step
- be wary of making the entry level simulation more complex (this is dynsim's primary advantage) - perhaps just make optional and show an example of how to do this?
- for Kalman filtering, we need to have read/write access to model states. How do we do this?
  - states are available as attributes of the classes; No point returning a vector of states from each object unless the order is specified/given. if documented properly (make sure all are listed), then could hand write methods of the model-class to get and put these by name. E.g. tank1.mLiquid .
  - or add a method "getStates" to ALL the component classes to be able to do this without figuring out which attributes are genuine states
  - write an EKF class and test these thoughts.
Tools
- simulation control dialog (start/pause/continue/stop, simulation time, dt?)
- simple user interface (table) for changing nominated variables during simulation (while sim stopped?)
- ability to read in a sequence of dictionaries or scripts and run simulations for each different scenario.
More block classes
- ????control.bin to simulate simple fast bin/tank level with no overflow (compare with 2-phase flow "process.tank" class)
  - or maybe just use integrate class like: intobj.run(dt,meas=(inflow-outflow),Ti=pct_per_tonne/60minperhour) ? (Ti calc needs checking - for example only, appropriately naming "intobj" as (say) "bin" would make code quite readable
- bin class which represents inventory/level as well as the plug flow of solids through the bin (of a vector of information such as grade or size)
- DCU building blocks?
- add derivative filtering on PID class
- arbitrary impulse response (FIR) ( vector of (time,value) tuples?)
- non-linear (hysteresis/stiction/rate-of-change-limit etc.; put in one block?, stiction can go in valve class as well)
- timer - two modes (pulseoutput argument): returns True exactly once after time is set, or returns True until timeout ( timerobj.run(dt,settime=None,pulseoutput=False) ). This will be useful for (amoung other things) putting together a simple finite state machine in a simulation, without having to create a whole FSM class.
- (Extended) Kalman Filter (for EKF, need to embed single timestep model into a class, ready for linearisation)
  - also consider UKF as an extra
- NN learn and execute class (using bpnn.py; note that scipy_bpnn.py, using numpy, is slower for small and medium sized networks)
- MPC algorithm (just start with unconstrained matrix solution)
  - how to solve without parameterising control actions
    - possibly use delta actions
      - these should be independent
      - independence makes easy
  - solve unconstrained case for initial values
    - just a matrix equation!
    - then use nonlinear minimisation
  - use full vector of PV/CO time values
    - more info to optimiser
  - use NN to learn control function (NLMPC)
    - fast execution
    - gen data by exhaustive search of space?
- process classes
  - pump (improved pump curve: flow vs pressure, speed, viscosity)
  - thickener model (mostly implemented, untested)
More example simulations (especially quirky ones!)
- constraint control (PID, model, or MPC based?)
- control of inverse response dynamics (build with two "lag"s; PID controller, DMC ...?)
create on-line demo simulations
- through a macro which collects form-based parameters, or collects parameters using a modified "PageComment" macro posting of a dictionary script (cut from elsewhere on the page, pasted and modified), or a special parser for an editable code block (logged in users only).
- runs the code (with a special flag to produce a .png file) and shows the plot on the page
- watch out for security! (ability to post and run arbitrary code)
develop a clean way to integrate SimPy classes to do hybrid simulation

Not To Do

? sequencing classes ? (flip-flops, state-machine tools etc - this may be better just done with if-then statements, i.e. too specialised to warrant dedicated classes)

CategoryActive