Background

Comer Duncan recently sent me an email asking how to translate models written for Berkeley Madonna into Modelica. He was specifically interested in some metabolism models developed by Kevin Hall at NIH.

Instead of simply answering this question via an email discussion, I thought it would be good to write about this issue. As is often the case, this kind of undertaking risks falling into the "Your Baby is Ugly" trap, although I hope that won't be the case here.

For this article, I am going to focus exclusively on a "straight translation" approach. This is a naive but insightful approach. In part two, I'll examine a better, but more involved, approach the brings to light many of the advantages that Modelica has to offer.

Straight Translation

These models contain several different types of statements. Let's take a look at some representative statements in Berkeley Madonna and then show what a literal translation to Modelica would look like.

System Parameters

The first type of statement are ones that set the value of so-called "system parameters" in Berkeley Madonna. Examples of system parameters include time step, tolerances, and start time. For example, to set the start time to 0, a Berkeley Madonna model would include the following line:

STARTTIME = 0

Other system parameters include STOPTIME, DTMIN, DTMAX, TOLERANCE and DTOUT.

In Modelica, the experimental conditions are associated with the model via a standard annotation. This helps segregate meta-data about how the problem is to be solved from the actual mathematical equations associated with the model. So in Modelica, we would indicate start time as follows:

  annotation(experiment(StartTime=0.0))

Differential Equations

Berkeley Madonna allows differential equations to be defined using the following syntax:

d/dt(R) = ...
R' = ...
FLOW R = ...

The idea here is that for some reservoir, R, the right hand side of these equations would define the rate at which the reservoir value changes with respect to time. The fact that such variables are called reservoirs shows the influence that Forrester's System Dynamics has on the approach in Berkeley Madonna. This isn't a criticism. Forrester's approach to system dynamics, with its concepts of stocks and flows, is quite intuitive. As we will see later, the formalism used in Modelica is a superset of Forrester's system dynamics (and several other formalisms).

In Modelica, a differential equation has only one representation:

der(R) = ...;

where der is a built-in operator. It is worth noting that Modelica is a declarative programming language, not an imperative programming language. As such, there is no "assignment" or directionality in a Modelica equation, only a relationship between two quantities. So the same equation could be represented in the following, completely equivalent, form:

... = der(R);

It is worth noting that Berkeley Madonna also allows the following forms, in order to be compatible with the STELLA modeling language:

R(t) = R(t-dt) + (...)*dt
R = R + dt*(...)

In the Berkeley Madonna User's Guide, it says these forms are not recommended "since the notation is more error-prone". At the risk of offending STELLA advocates, I'm afraid I have to agree. It isn't just the potential complexity of such expressions, it is the fact that it implicitly imposes a (forward Euler) solution method on the model which is completely unnecessary (and not a particularly good choice at that). Mixing the "problem statement" and the "solution method" is avoided as much as possible in the Modelica approach so there are no equivalent forms in Modelica. That is not to say that discrete equations (e.g., z-transforms) are not allowed in Modelica. But such equations must be represented in terms of some underlying "clock" and not in terms of solver time steps.

Berkeley Madonna allows for higher-derivatives to be expressed in the language through the use of multiple "primes", e.g., u', u''. Modelica does not allows this. Differential equations are restricted to first-order and intermediate variables must be introduced to represent higher-order derivatives.

Initialization

For each reservoir, it is necessary to specify an initial condition. In Berkeley Madonna, this can be done in the following ways:

INIT R = ...
INIT (R) = ...
R(STARTIME) = ...

where the right hand side represents the value the reservoir should have at the start of the simulation.

Initialization in Modelica is actually quite a rich topic. For since we are only concerned with straight translation at the moment, the equivalent in Modelica would be:

initial equation
  R = ...;

Comments

The comment syntax in Berkeley Madonna has two forms. The first form is a "curly-bracket comment" which can span multiple lines, e.g.,

{This is a comment}
{This is
 also a comment}

The second form is a single-line comment which starts with a ; and is terminated by the end of the line, e.g.,

; This is a comment

Modelica has both of these types of comments with slightly different syntax:

/* This is
   a multi-line
   comment */
// This is a single-line comment

But in addition, Modelica has "descriptive strings" that can be associated with entities in the language. These are not comments because they are semantically associated with different elements. Such descriptive strings can then be used in parameter dialogs and other user interface elements. For example, if I declare a variable as follows:

parameter Real Na_b=4000 "Baseline sodium intake in mg/d";

the string "baseline sodium intake in mg/d" is not a comment but rather a description of the variable Na_b and this description can be used when interacting or documenting the model.

Furthermore, Modelica allows you to formally associated physical units with a quantity. As such, an even better representation in Modelica would be:

type DailyRate = Real(final unit="mg/d");
parameter DailyRate Na_b=4000 "Baseline sodium intake";

In this way, we can leave the units out of the description because the variable is declared to be of type DailyRate, which already indicates that the units are in milligrams per day. Even better, the Modelica specification defines a grammar for these unit definitions so that the consistency of units in equations can be checked. This goes well beyond just commenting or documenting the model because it allows unit errors to be automatically detected by tools which can catch a lot potential errors!

Variables and Equations

In Berkeley Madonna, models consist of a list of equations. As we have seen already, some of these equations are for "system parameters". Others are for variables in our models. Examples of equations include:

ECP = 0.732*BM + 0.01087*ECW_b
d/dt (Lipol_diet) = (Lipol_diet_target - Lipol_diet)/tau_lip
init Lipol_diet = 1
Kurine = IF Ketogen < Kspill THEN 0 ELSE Ketogen*KUmax/(KGmax-Kspill)-KUmax/(KGmax/Kspill-1)

In these cases, you have either a variable name or the derivative of a variable on the left hand side and some kind of expression on the right hand side. With only minor syntactic differences, we can represent these same equations in Modelica:

equation
  ECP = 0.732*BM + 0.01087*ECW_b;
  der(Lipol_diet) = (Lipol_diet_target - Lipol_diet)/tau_lip;
  Kurine = if Ketogen < Kspill then 0 else Ketogen*KUmax/(KGmax-Kspill)-KUmax/(KGmax/Kspill-1);
initial equation
  Lipol_diet = 1;

The main differences are the termination of equations with semicolons, the use of the der operator instead of d/dt and the fact that Modelica is case sensitive (if vs IF).

However, there is another important difference with regard to variables that we mentioned already in the section on commenting, which is that variables in Modelica must be declared. So a more complete fragment of Modelica code for the previous three equations would be:

model MetabolismModel
  Real ECP "Extracellular protein";
  Real Lipol_diet;
  Real Kurine;
  ...
equation
  ECP = 0.732*BM + 0.01087*ECW_b "Wang AJCN 2003";
  der(Lipol_diet) = (Lipol_diet_target - Lipol_diet)/tau_lip;
  if Ketogen < Kspill then
    Kurine = 0;
  else
    Kurine = Ketogen*KUmax/(KGmax-Kspill)-KUmax/(KGmax/Kspill-1);
  end if;
  ...
initial equation
  Lipol_diet = 1;
  ...
end MetabolismModel;

Here we see the mostly complete text of a model. We see both the declarations of the variables (indicating the type of the variable along with an optional description) as well as the equations associated with the variable.

Isn't the Berkeley Madonna syntax simpler? Perhaps. For simple problems it might seem like an advantage to have such a simple syntax. But for complex problems, using some explicit syntax to help convey your overall intent can go a long way toward providing better diagnostic error messages and catching errors.

At the risk of jumping ahead a little bit, it is worth pointing out that building complex models in Modelica would not be done this way (i.e. declaring lots of variables and lots of equations all in a "flat" file like this). But for now, we will remain focused on a straight translation.

Datasets/Tables

One interesting thing that Berkeley Madonna has is the notion of datasets. These datasets are really a combination of two things. The first is the underlying data (presumably represented on on some multi-dimensional regular grid, although I didn't confirm that). The other is a bunch of implicitly defined functions for interpolating over the data in the dataset. These implicit functions are identified by the name of the dataset preceded by a # character, e.g., #temperature(...), #R20BW(...).

Modelica (or more specifically, the Modelica Standard Library) includes a collection of table models that are similar, but not exactly equivalent.

An important caveat here is that in practice it is often the case that you might wish to use a table in some cases, a set of mathematical expressions in another and perhaps even a nested sub-model with its own states in another. We'll talk later about how Modelica can accommodate all of these uses in a framework that is "type safe".

Automatic Translation

It turns out, it is pretty straight forward to translate Kevin Hall's models (which have only appeared in fragments in this article) into Modelica code. To do this, I wrote a relatively simply Python script although I should point out that if I needed something that was "production quality", I would use a real lexer and parser. It would take me perhaps twice as long to write, but it would be infinitely more reliable and robust.

The following Python code is almost sufficient to translate the model by Kevin Hall that I mentioned at the beginning into Modelica¹:

#!/usr/bin/python
import re

def _process_eq(groups, eqs, line):
    lhs = groups[0]
    rhs = groups[1].split(";")
    if len(rhs)==1:
        eqs.append({"var": lhs, "expr": rhs[0].strip(),
                    "desc": None})
    elif len(rhs)>1:
        eqs.append({"var": lhs, "expr": rhs[0].strip(),
                    "desc": ";".join(rhs[1:])})
    else:
        print "Unable to parse equation: ", line

class Translator:
    EMPTY = re.compile("^\s*(;.*)?$")
    ASSIGNMENT = re.compile("\s*(\w+)\s*=\s*(.*)")
    INIT = re.compile("\s*[Ii][Nn][Ii][Tt]\s*(\w+)\s*=(.*)")
    DIFFEQ = re.compile("\s*d\/dt\s*\(?(\w+)\)?\s*=\s*(.*)")
    CCOMMENT = re.compile("\s*{([^}]*)}\s*")
    def __init__(self, name, f):
        self.name = name
        self.equations = []
        self.diffequations = []
        self.initequations = []
        self.file = f
        self._parse()
    def _parse(self):
        lines = self.file.readlines()
        for line in lines:
            self._parseline(line)
    def _parseline(self, line):
        if self.EMPTY.match(line):
            return
        if self.CCOMMENT.match(line):
            return
        if self.DIFFEQ.match(line):
            match = self.DIFFEQ.match(line)
            _process_eq(match.groups(), self.diffequations, line)
        elif self.INIT.match(line):
            match = self.INIT.match(line)
            _process_eq(match.groups(), self.initequations, line)
        elif self.ASSIGNMENT.match(line):
            match = self.ASSIGNMENT.match(line)
            _process_eq(match.groups(), self.equations, line)
        else:
            print "?", line
        return
    def render(self):
        experiment = {}
        # Extract experimental settings
        equations = list(self.equations)
        for eq in equations:
            var = eq["var"]
            exp = eq["expr"]
            if var=="STARTTIME":
                experiment["StartTime"] = exp;
                self.equations.remove(eq)
            if var=="STOPTIME":
                experiment["StopTime"] = exp;
                self.equations.remove(eq)
            if var=="TOLERANCE":
                experiment["Tolerance"] = exp;
                self.equations.remove(eq)
            if var=="DTMIN" or var=="DTMAX" or var=="DTOUT":
                self.equations.remove(eq)
        print "model %s" % (self.name)
        for eq in self.equations:
            if hasattr(eq,"desc"):
                print """  Real %s "%s";""" % (eq["var"], eq["desc"])
            else:
                print """  Real %s;""" % (eq["var"])
        print "initial equation"
        for ieq in self.initequations:
            print """  %s=%s;""" % (ieq["var"], ieq["expr"])
        print "equation"
        for eq in self.equations:
            print """  %s=%s;""" % (eq["var"], eq["expr"])
        for deq in self.diffequations:
            print """  der(%s)=%s;""" % (deq["var"], deq["expr"])
        print "end %s;" % (self.name)

fp = open("hallcode.m", "r")
t = Translator("HallModel", fp)
t.render()

For the curious, I have included the resulting Modelica model (which required just a handful of manual tweaks).

Now What?

I've mentioned a few times that this kind of straight translation is really not the way to go about this. So far, what I've shown is a bit naive, but it does work. In part two of this article, we'll look at how you could improve on this approach.